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