当前位置：首页 > news >正文

C#，数值计算，求解微分方程的预测校正法（修正欧拉法）算法与源代码

news 来源：原创 2024/5/20 11:39:16

Leonhard Euler

1 微分方程

微分方程，是指含有未知函数及其导数的关系式。解微分方程就是找出未知函数。
微分方程是伴随着微积分学一起发展起来的。微积分学的奠基人Newton和Leibniz的著作中都处理过与微分方程有关的问题。微分方程的应用十分广泛，可以解决许多与导数有关的问题。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。此外，微分方程在化学、工程学、经济学和人口统计等领域都有应用。
数学领域对微分方程的研究着重在几个不同的面向，但大多数都是关心微分方程的解。只有少数简单的微分方程可以求得解析解。不过即使没有找到其解析解，仍然可以确认其解的部分性质。在无法求得解析解时，可以利用数值分析的方式，利用电脑来找到其数值解。动力系统理论强调对于微分方程系统的量化分析，而许多数值方法可以计算微分方程的数值解，且有一定的准确度。

2 预估-校正法

对于给定的带有初始条件y(x0)=y0的微分方程dy/dx=f(x, y)，用预测-校正法求近似解。
预估-校正法: 预估-校正法又称修正-欧拉法。
在欧拉方法中，在一个点绘制切线，并计算给定步长的斜率。
因此，这种方法对于线性函数效果最好，但是对于其他情况，仍然存在截断误差。
为了解决这个问题，引入了修正欧拉法。

3 源程序

using System;
using System.Collections;
using System.Collections.Generic;

namespace Legalsoft.Truffer.Algorithm
{
/// <summary>
/// 给定微分方程的一阶偏导方程
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
/// <returns></returns>
public delegate double SDE_Equation(double x, double y);

/// <summary>
/// 求解微分方程的预测校正法或修正欧拉法
/// </summary>
public static partial class Algorithm_Gallery
{
public static SDE_Equation dydx = null;

private static double Predict(double x, double y, double h)
{
double y1p = y + h * dydx(x, y);
return y1p;
}

private static double Correct(double x, double y, double x1, double y1, double h)
{
double e = 0.00001;
double y1c = y1;

do
{
y1 = y1c;
y1c = y + 0.5 * h * (dydx(x, y) + dydx(x1, y1));
}
while (Math.Abs(y1c - y1) > e);

return y1c;
}

public static double DE_Modified_Euler_Solve(double x, double xn, double y, double h)
{
while (x < xn)
{
double x1 = x + h;
double y1p = Predict(x, y, h);
double y1c = Correct(x, y, x1, y1p, h);
x = x1;
y = y1c;
}
return y;
}
}
}

POWER BY TRUFFER.CN

4 源代码

using System;
using System.Collections;
using System.Collections.Generic;namespace Legalsoft.Truffer.Algorithm
{/// <summary>/// 给定微分方程的一阶偏导方程/// </summary>/// <param name="x"></param>/// <param name="y"></param>/// <returns></returns>public delegate double SDE_Equation(double x, double y);/// <summary>/// 求解微分方程的预测校正法或修正欧拉法/// </summary>public static partial class Algorithm_Gallery{public static SDE_Equation dydx = null;private static double Predict(double x, double y, double h){double y1p = y + h * dydx(x, y);return y1p;}private static double Correct(double x, double y, double x1, double y1, double h){double e = 0.00001;double y1c = y1;do{y1 = y1c;y1c = y + 0.5 * h * (dydx(x, y) + dydx(x1, y1));}while (Math.Abs(y1c - y1) > e);return y1c;}public static double DE_Modified_Euler_Solve(double x, double xn, double y, double h){while (x < xn){double x1 = x + h;double y1p = Predict(x, y, h);double y1c = Correct(x, y, x1, y1p, h);x = x1;y = y1c;}return y;}}
}