实验八 一阶常微分方程初值问题Matlab实现

不实名感谢某位大兄弟的贡献

【作业内容】: 用Euler法、改进的欧拉法、四阶龙格-库塔法求如下一阶常微分方程

{y′=y−2xyy(0)=1\left\{\begin{array}{l} {y}^{\prime}={y}-\frac{2 {x}}{{y}} \\ {y}({0})={1} \end{array}\right. {y′=y−y2x​y(0)=1​

x∈[0,2]\mathbf{x} \in[0,2] x∈[0,2]

初值问题的解（小数点后保留8位小数），取h = 0.05.

【作业要求】：

将各种方法求出的解连同精确解使用同一个表格进行数据比较；画出精确解（见教材）的图形，并且将其它3种方法得到的近似解画出图形,进行比较，见如上的例子7.1。解释各种求解方法。

【题目分析及所采用的算法】

按题目要求，需要我们计算一阶常微分方程问题，分别使用Euler法、改进的欧拉法、四阶龙格-库塔法来进行迭代计算。

Euler法：

{yi+1=yi+hf(xi,yi)y0=y(x0),i=0,1,…,n\left\{\begin{array}{l}{y}_{\mathrm{i}+1}={y}_{\mathrm{i}}+{h} {f}\left({x}_{\mathrm{i}}, {y}_{\mathrm{i}}\right) \\ {y}_{0}={y}\left({x}_{0}\right)\end{array}, {i}={0}, {1}, \ldots, {n}\right. {yi+1​=yi​+hf(xi​,yi​)y0​=y(x0​)​,i=0,1,…,n

hhh为步长。依次从初值y0=y(x0)y_0 = y(x_0)y0​=y(x0​)开始迭代。

改进的欧拉法(将梯形公式和Euler法结合，精度更高):

yp=yi+hf(xi,yi)yc=yi+hf(xi+1,yp)yi+1=12(yp+yc)i=0,1,2…,n−1\begin{array}{l}{y}_{{p}}={y}_{{i}}+{h f}\left({x}_{{i}}, {y}_{{i}}\right) \\ {y}_{{c}}={y}_{{i}}+{h f}\left({x}_{{i}+{1}}, {y}_{{p}}\right) \\ {y}_{{i}+{1}}=\frac{{1}}{{2}}\left({y}_{{p}}+{y}_{{c}}\right) \\ {i}={0}, {1}, {2} \ldots, {n}-{1}\end{array} yp​=yi​+hf(xi​,yi​)yc​=yi​+hf(xi+1​,yp​)yi+1​=21​(yp​+yc​)i=0,1,2…,n−1​

hhh为步长。依次从初值y0=y(x0)y_0 = y(x_0)y0​=y(x0​)开始迭代。

四阶龙格-库塔法:

{k1=f(xi,yi)k2=f(xi+h2,yi+h2k1)k3=f(xi+h2,yi+h2k2)k4=f(xi+h,yi+hk3)yi+1=yi+h6(k1+2k2+2k3+k4)i=0,1,…,n\left\{\begin{array}{l}{k}_{1}={f}\left({x}_{i}, {y}_{i}\right) \\ {k}_{2}={f}\left({x}_{i}+\frac{{h}}{2}, {y}_{i}+\frac{{h}}{2} {k}_{1}\right) \\ {k}_{3}={f}\left({x}_{i}+\frac{{h}}{2}, {y}_{i}+\frac{{h}}{2} {k}_{2}\right) \\ {k}_{4}={f}\left({x}_{i}+{h}, {y}_{i}+{h} {k}_{3}\right) \\ {y}_{{i}+{1}}={y}_{{i}}+\frac{{h}}{{6}}\left({k}_{{1}}+{2} {k}_{{2}}+{2} {k}_{3}+{k}_{4}\right) \quad {i}={0}, {1}, \ldots, {n}\end{array}\right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​k1​=f(xi​,yi​)k2​=f(xi​+2h​,yi​+2h​k1​)k3​=f(xi​+2h​,yi​+2h​k2​)k4​=f(xi​+h,yi​+hk3​)yi+1​=yi​+6h​(k1​+2k2​+2k3​+k4​)i=0,1,…,n​

hhh为步长。依次从初值y0=y(x0)y_0 = y(x_0)y0​=y(x0​)开始迭代。

【主要代码（要求程序代码有详细的注释）】

代码含有相关注释。

%Ordinary_differential.mformat long;% 步长h = 1/20;% 右边界xf = 2;% 设置步长x = 0 : h : xf;% 已知函数Func = @(x, y)(y - ((2 * x) / y));%Euler% 初始化y数组y_E = zeros(1, length(x));y_E(1) = 1;% 复现公式 & 迭代计算for i = 1:(length(x)-1)y_E(i + 1) = y_E(i) + Func(x(i), y_E(i)) * h;end%Improve Euler% 初始化y数组y_IE = zeros(1, length(x));y_IE(1) = 1;% 复现公式 & 迭代计算for i = 1:(length(x) - 1)p = y_IE(i) + h * Func(x(i), y_IE(i));c = y_IE(i) + h * Func(x(i + 1), p);y_IE(i + 1) = (1/2) * (p + c);end%Rung Kutta% 初始化y数组y_RK = zeros(1, length(x));y_RK(1) = 1;% 复现公式 & 迭代计算for i=1:(length(x) - 1) k1 = Func(x(i), y_RK(i));k2 = Func(x(i) + 0.5 * h, y_RK(i) + 0.5 * h * k1);k3 = Func((x(i) + 0.5 * h), (y_RK(i) + 0.5 * h * k2));k4 = Func((x(i) + h), (y_RK(i) + k3 * h));y_RK(i + 1) = y_RK(i) + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4) * h;end%Exact% 使用ode45计算精确函数x0 = 0;y0 = 1;yspan = [x0 xf];[x_ode45, y_ode45] = ode45(Func,yspan,y0);% 绘图subplot(221)plot(x_ode45, y_ode45, '-'); xlabel('x');ylabel('y');legend('Exact');subplot(222)plot(x, y_E, '-');xlabel('x');ylabel('y');legend('Euler');subplot(223)plot(x, y_IE, '-');xlabel('x');ylabel('y');legend('Improve Euler');subplot(224)plot(x, y_RK, '-');xlabel('x');ylabel('y');legend('Rung Kutta');% tableres = [round(transpose(x), 8), round(transpose(y_E), 8), round(transpose(y_IE), 8), round(transpose(y_RK), 8), round(y_ode45, 8)];table = array2table(res, 'VariableNames', {'x','Euler','Improve Euler', 'Rung Kutta', 'Exact'});table;

【结果比较与（或）分析】

绘图：

明显的看出来欧拉法并没有其他两种方法得来的曲线光滑，更不精确。

可以从数值上进一步判断其精确程度：

>> tabletable =41×5 tablex Euler Improve Euler Rung KuttaExact ____ __________ _____________ __________ __________0 1 11 10.051.051.048869051.04880886 1.048810180.11.09773811.095561121.09544514 1.095445380.15 1.143515361.140344591.14017546 1.140174580.2 1.187573681.18343698 1.183216 1.183216060.25 1.230111311.225017491.22474493 1.224745260.3 1.271293511.265235851.26491113 1.264911110.35 1.311260171.304218771.30384056 1.303840290.41.35013131.342074541.34164088 1.341640950.45 1.388011111.37889671.37840498 1.378405090.5 1.424991181.414766631.41421368 1.414213620.551.46115281.449755741.44913781 1.449137690.6 1.496568931.483927111.48323985 1.483239910.65 1.531305671.517336791.51657526 1.516575290.7 1.565423521.550034921.54919353 1.549193450.75 1.598978361.582066531.58113904 1.581138970.8 1.63331.613472331.61245178 1.612451820.85 1.664604511.644289251.64316793 1.643167930.9 1.696771561.674550961.67332034 1.673320260.95 1.728568231.70428831.70293895 1.702938891 1.760037861.733529621.73205115 1.732051171.05 1.791222791.762301091.76068206 1.760682041.1 1.822164761.790626971.78885479 1.788854711.15 1.852905241.818529811.816590661.81659061.21.88348581.846030711.84390938 1.843909391.25 1.913948441.873149421.87082923 1.870829191.3 1.944335841.899904561.89736719 1.897367091.35 1.974691761.926313731.92353905 1.923538981.4 2.005061251.952393631.94935958 1.949359561.45 2.035491011.97816021.97484254 1.974842471.5 2.066029672.00362872.00000085 2.000000731.55 2.096728132.0288138 2.02484662.02484651.6 2.127639842.053729732.04939117 2.049391111.65 2.158821132.078390272.07364525 2.073645151.7 2.190331592.102808922.09761891 2.097618771.75 2.222234352.126998932.12132168 2.121321541.82.25459652.150973392.14476252 2.144762411.85 2.287489422.174745292.16794994 2.167949781.92.32098922.198327642.19089198 2.190891781.95 2.355177012.221733482.21359628 2.213596082 2.390139542.244976012.236070082.2360699>>

可以看出来。精确程度：欧拉法 < 改进的欧拉法 < 龙格库塔法 。