实验目的及要求：

1、熟悉各种初值问题的算法，编出算法程序；

2、明确各种算法的精度寓所选步长有密切关系；

3、熟悉在Matlab平台上直接求解常微分方程初值问题。

实验内容：

1、编写改进的欧拉公式通用子程序，取

求解初值问题 。

2、编写四阶经典龙格-库塔方法子程序，取

求解初值问题，并计算y(0.4)的近似值，要求小数点后保留4位有效数字。

实验步骤与程序：

1. 改进的欧拉格式：

1. 算法分析：

1. 改进的欧拉格式流程图：

或者

1. 改进欧拉格式的MATLAB主程序

被调用的ou_la_gaijin.m文件

function [x,y]=ou_la_gaijin(x0,y0,h,N,ydot_fun)

x=zeros(1,N+1);

y=zeros(1,N+1);

x(1)=x0;

y(1)=y0;

for n=1:N

    x(n+1)=x(n)+h;

    ybar=y(n)+h*feval(ydot_fun,x(n),y(n));

    y(n+1)=y(n)+(h/2)*(feval(ydot_fun,x(n),y(n))+feval(ydot_fun,x(n+1),ybar));

运行的ou_la_gaijin1.m文件

ydot_fun=inline('1/(1+x^2)-2*y^2');

[x,y]=ou_la_gaijin(0,0,0.1,40,ydot_fun)

1. 经典四阶龙格-库塔格式：

1. 算法分析

1. 经典四阶龙格-库塔格式流程图：

1. 四阶龙格-库塔格式的MATLAB主程序

被调用的long_ge_ku_ta_4j.m文件

function [x,y]=long_ge_ku_ta_4j(x0,y0,h,N,fun)

x=zeros(1,N+1);

!生成矩阵

y=zeros(1,N+1);

x(1)=x0;

y(1)=y0;

for n=1:N

    x(n+1)=x(n)+h;

    k1=feval(fun,x(n),y(n));

    k2=feval(fun,x(n)+h/2,y(n)+(h*k1)/2);

    k3=feval(fun,x(n)+h/2,y(n)+(h*k2)/2);

    k4=feval(fun,x(n)+h,y(n)+h*k3)

    y(n+1)=y(n)+(h*(k1+2*k2+2*k3+k4))/6;

end

fun=@(x,y)8-3*y;

 [x,y]=long_ge_ku_ta_4j(0,2,0.2,5,fun)

运行的long_ge_ku_ta_4j_1.m文件

fun=@(x,y)8-3*y;

!@()的意思是定义函数

[x,y]=long_ge_ku_ta_4j(0,2,0.2,5,fun)

运行结果：

结果分析与讨论：

1. 改进欧拉法是对欧拉算法的改进方法。微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值，这个过程称为离散化。实现离散化的基本途径是用向前差商来近似代替导数，这就是欧拉算法实现的依据。欧拉(Euler)算法是数值求解中最基本、最简单的方法，但其求解精度较低，
2. 改进欧拉法是先用欧拉法求得一个初步的近似值，称为预报值，然后用它替代梯形法右端的yi+1再直接计算fi+1，得到校正值yi+1，这样建立的预报-校正系统称为改进的欧拉格式：预报值 y~i+1=yi + h*f(xi,yi)  校正值 yi+1 =yi+(h/2)*[f(xi,yi)+f(xi+1,y~i+1)]它有下列平均化形式：yp=yi+h*f(xi,yi)且 yc=yi+h*f(xi+1,yp)且 yi+1=(yp+yc)/2，它的局部截断误差为O(h^3)，可见，改进欧拉格式较欧拉格式提高了精度，其截断误差比欧拉格式提高了一阶。
3.  四阶龙格-库塔法算法中，k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。