实现了一些常见的带初值的常微分方程的算法。
/// <summary> /// 欧拉迭代公式,求出区间(x0,x]中每隔步长h的精度为e的近似值 /// </summary> /// <param name="fun">fun为x的函数即 dy/dx=fun(x,y)</param> /// <param name="N">将区间分为N段</param> /// <param name="x0">f(x0)=y0</param> /// <param name="y0">f(x0)=y0</param> /// <param name="x">所求区间的右端点</param> /// <param name="e">迭代精度</param> /// <returns>返回区间每隔h的函数值</returns> public static double[] Euler_e(Func<double, double, double> fun, int N, double x0, double y0, double x, double e = 1e-10) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h, yn); double res_e = 0; while (Math.Abs(res_e - res[i]) >= e) { res_e = res[i]; res[i] = yn + h / 2 * (fun(x0 - h, yn) + fun(x0, res_e)); } } return res; } /// <summary> /// 欧拉迭代公式,求出区间(x0,x]中每隔步长h的精度为e的近似值 /// </summary> /// <param name="fun">fun为x的函数即 dy/dx=fun(x)</param> /// <param name="N">将区间分为N段</param> /// <param name="x0">f(x0)=y0</param> /// <param name="y0">f(x0)=y0</param> /// <param name="x">所求区间的右端点</param> /// <param name="e">迭代精度</param> /// <returns>返回区间每隔h的函数值</returns> public static double[] Euler_e(Func<double, double> fun, int N, double x0, double y0, double x, double e = 1e-10) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h); double res_e = 0; while (Math.Abs(res_e - res[i]) >= e) { res_e = res[i]; res[i] = yn + h / 2 * (fun(x0 - h) + fun(x0)); } } return res; } //欧拉预估-矫正公式 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x,y) //f(x0)=y0 public static double[] Euler_pre(Func<double, double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h, yn); res[i] = yn + h / 2 * (fun(x0 - h, yn) + fun(x0, res[i])); } return res; } //欧拉预估-矫正公式 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x) //f(x0)=y0 public static double[] Euler_pre(Func<double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h); res[i] = yn + h / 2 * (fun(x0 - h) + fun(x0)); } return res; } //二阶龙格-库塔算法 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x,y) //f(x0)=y0 public static double[] R_K2(Func<double, double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0, yn); double k2 = h * fun(x0 + 2.0 / 3 * h, yn + 2.0 / 3 * k1); res[i] = yn + 1.0 / 4 * (k1 + 3 * k2); x0 += h; } return res; } //二阶龙格-库塔算法 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x) //f(x0)=y0 public static double[] R_K2(Func<double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0); double k2 = h * fun(x0 + 2.0 / 3 * h); res[i] = yn + 1.0 / 4 * (k1 + 3 * k2); x0 += h; } return res; } //四阶龙格-库塔算法 public static double[] R_K4(Func<double, double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0, yn); double k2 = h * fun(x0 + 0.5 * h, yn + 0.5 * k1); double k3 = h * fun(x0 + 0.5 * h, yn + 0.5 * k2); double k4 = h * fun(x0 + h, yn + k3); res[i] = yn + 1.0 / 6 * (k1 + 2 * (k2 + k3) + k4); x0 += h; } return res; } //四阶龙格-库塔算法re public static double[] R_K4(Func<double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0); double k2 = h * fun(x0 + 0.5 * h); double k3 = h * fun(x0 + 0.5 * h); double k4 = h * fun(x0 + h); res[i] = yn + 1.0 / 6 * (k1 + 2 * (k2 + k3) + k4); x0 += h; } return res; }