C#，数值计算，求解微分方程的吉尔（Gear）四阶方法与源代码
1 微分方程
微分方程，是指含有未知函数及其导数的关系式。解微分方程就是找出未知函数。
微分方程是伴随着微积分学一起发展起来的。微积分学的奠基人Newton和Leibniz的著作中都处理过与微分方程有关的问题。微分方程的应用十分广泛，可以解决许多与导数有关的问题。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。此外，微分方程在化学、工程学、经济学和人口统计等领域都有应用。
数学领域对微分方程的研究着重在几个不同的面向，但大多数都是关心微分方程的解。只有少数简单的微分方程可以求得解析解。不过即使没有找到其解析解，仍然可以确认其解的部分性质。在无法求得解析解时，可以利用数值分析的方式，利用电脑来找到其数值解。 动力系统理论强调对于微分方程系统的量化分析，而许多数值方法可以计算微分方程的数值解，且有一定的准确度。
2 数值解法
作为数值分析的基础内容，常微分方程数值解法的研究已发展得相当成熟，理论上也颇为完善，各类有实用价值的算法已经建立，并已形成计算机软件。它处理问题的思路与方法常可用于偏微分方程的数值求解。主要研究以下三类定解问题的数值解法：初值问题、两点边值问题与特征值问题。初值问题的数值解法应用广泛，是常微分方程数值解法的主要内容。在这方面有突出贡献的学者当推达赫奎斯特（Dahlquist，G.）、巴特赫尔（Butcher，J.C.）及吉尔（Gear，C.W.）等人。两点边值问题及特征值问题的研究相对较为薄弱，其中凯勒尔（Keller，H.B.）的工作影响较大。
Gear, C.William
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>
/// 求解微分方程的吉尔四阶方法
/// C# program to implement Gill's method
/// </summary>
public static partial class Algorithm_Gallery
{
public static SDE_Equation dydx = null;
/// <summary>
/// 求解微分方程的吉尔四阶方法
/// </summary>
/// <param name="x0">起点x坐标</param>
/// <param name="y0">起点y坐标</param>
/// <param name="x">求值点x坐标</param>
/// <param name="step">步长</param>
/// <returns></returns>
public static double SDE_Gill_Method(double x0, double y0, double x, double step)
{
int n = (int)((x - x0) / step);
double y = y0;
for (int i = 1; i <= n; i++)
{
double k1 = step * dydx(x0, y);
double k2 = step * dydx(x0 + 0.5 * step, y + 0.5 * k1);
double k3 = step * dydx(x0 + 0.5 * step, y + 0.5 * (-1 + Math.Sqrt(2)) * k1 + k2 * (1 - 0.5 * Math.Sqrt(2)));
double k4 = step * dydx(x0 + step, y - (0.5 * Math.Sqrt(2)) * k2 + k3 * (1 + 0.5 * Math.Sqrt(2)));
y = y + (1.0 / 6) * (k1 + (2 - Math.Sqrt(2)) * k2 + (2 + Math.Sqrt(2)) * k3 + k4);
x0 = x0 + step;
}
return y;
}
}
}
使用该方法的参考代码（POWER BY 315SOFT.COM）：
using Legalsoft.Truffer.Algorithm;
namespace Legalsoft.Drive
{
public partial class Form1 : Form
{
public double func(double x, double y)
{
return x/2 + y*y;
}
private void button1_Click(object sender, EventArgs e)
{
Algorithm_Gallery.dydx = func;
MessageBox.Show("result="+ Algorithm_Gallery.SDE_Gill_Method(0.0,0.0,0.5,30));
}
}
}
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>/// 求解微分方程的吉尔四阶方法/// C# program to implement Gill's method/// </summary>public static partial class Algorithm_Gallery{public static SDE_Equation dydx = null;/// <summary>/// 求解微分方程的吉尔四阶方法/// </summary>/// <param name="x0">起点x坐标</param>/// <param name="y0">起点y坐标</param>/// <param name="x">求值点x坐标</param>/// <param name="step">步长</param>/// <returns></returns>public static double SDE_Gill_Method(double x0, double y0, double x, double step){int n = (int)((x - x0) / step);double y = y0;for (int i = 1; i <= n; i++){double k1 = step * dydx(x0, y);double k2 = step * dydx(x0 + 0.5 * step, y + 0.5 * k1);double k3 = step * dydx(x0 + 0.5 * step, y + 0.5 * (-1 + Math.Sqrt(2)) * k1 + k2 * (1 - 0.5 * Math.Sqrt(2)));double k4 = step * dydx(x0 + step, y - (0.5 * Math.Sqrt(2)) * k2 + k3 * (1 + 0.5 * Math.Sqrt(2)));y = y + (1.0 / 6) * (k1 + (2 - Math.Sqrt(2)) * k2 + (2 + Math.Sqrt(2)) * k3 + k4);x0 = x0 + step;}return y;}}
}