计算方法（四）带初值的常微分方程解法

 /// <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;
        }