改进欧拉方程求初值（十四）

改进欧拉方程的公式如下：

$y_p=y_n+hf\left(x_n,y_n\right)$
$y_c=y_n+hf\left(x_{n+1},y_n\right)$
$y_{n+1}=\frac{1}{2}\left(y_p+y_c\right)$

利用改进后的欧拉方程建立M函数文件求解初值。函数文件如下所示：

function [x,y] = adeuler(f,x0,y0,xf,h)
n = fix((xf-x0)/h);
y(1)=y0;
x(1)=x0;
for i=1:n
    x(i+1) = x0+i*h;
    yp = y(i) + h*feval(f,x(i),y(i));
    yc = y(i) + h*feval(f,x(i+1),yp);
    y(i+1) = (yp+yc)/2;%欧拉公式
end

例 求解初值条件$\{\begin{array}{l}{y}^{\prime }=y-\frac{2x}{y}(0<x<1)\\ y(0)=1\end{array}$

MATLAB的程序如下：

 [x,y] = adeuler('f',0,1,1,0.1)

x =

  1 至 10 列

​     0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000

  11 列

1.0000


y =

  1 至 10 列

1.0000    1.0959    1.1841    1.2662    1.3434    1.4164    1.4860    1.5525    1.6165    1.6782

  11 列

1.7379

y1 = (1+2*x).^0.5

y1 =

  1 至 10 列

1.0000    1.0954    1.1832    1.2649    1.3416    1.4142    1.4832    1.5492    1.6125    1.6733

  11 列

1.7321

plot(x,y,x,y1,'--')

图像结果如下所知：

由上面绘制曲线可以看出，改进的欧拉方程解法要比原来的方法优秀，数值解和解析解曲线能够重合。