实验四常微分方程初值问题数值解法讲解.docx
《实验四常微分方程初值问题数值解法讲解.docx》由会员分享，可在线阅读，更多相关《实验四常微分方程初值问题数值解法讲解.docx（13页珍藏版）》请在冰豆网上搜索。
实验四常微分方程初值问题数值解法讲解
实验报告
课程名称数值分析
实验项目常微分方程问题初值问题数值解法
一.实验目的
1、理解如何在计算机上实现用Euler法、改进Euler法、Runge－Kutta算法求一阶常微分方程初值问题
的数值解。
2、利用图形直观分析近似解和准确解之间的误差。
3、学会Matlab提供的ode45函数求解微分方程初值问题。
二、实验要求
（1）按照题目要求完成实验内容；
（2）写出相应的Matlab程序；
（3）给出实验结果（可以用表格展示实验结果）；
（4）分析和讨论实验结果并提出可能的优化实验。
（5）写出实验报告。
三、实验步骤
1、用编好的Euler法、改进Euler法计算书本P167的例1、P171例题3。
（1）取，求解初值问题
（2）取，求解初值问题
2、用Runge－Kutta算法计算P178例题、P285实验任务
（2）
（1）取，求解初值问题
（2）求初值问题
的解在处的近似值，并与问题的解析解相比较。
3、用Matlab绘图函数plot（x,y）绘制P285实验任务
（2）的精确解和近似解的图形。
4、使用matlab中的ode45求解P285实验任务
（2），并绘图。
四、实验结果
1、Euler算法程序、改进Euler算法程序；
Euler算法程序：
function[x,y]=euler_f（ydot_fun,x0,y0,h,N）
%Euler（向前）公式，其中
%ydot_fun---一阶微分方程的函数
%x0,y0---初始条件
%h---区间步长
%N---区间的个数
%x---Xn构成的向量
%y---Yn构成的向量
x=zeros（1,N+1）;y=zeros（1,N+1）;x
（1）=x0;y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
y（n+1）=y（n）+h*feval（ydot_fun,x（n）,y（n））;
end
改进Euler算法程序：
function[x,y]=euler_r（ydot_fun,x0,y0,h,N）
%改进Euler公式，其中
%ydot_fun---一阶微分方程的函数
%x0,y0---初始条件
%h---区间步长
%N---区间的个数
%x---Xn构成的向量
%y---Yn构成的向量
x=zeros（1,N+1）;y=zeros（1,N+1）;x
（1）=x0;y
（1）=y0;
forn=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））;
end
2、用Euler算法程序、改进Euler算法求解P167例题1的运行结果；
（1.）Euler算法程序：
x=
Columns1through8
00.10000.20000.30000.40000.50000.60000.7000
Columns9through11
0.80000.90001.0000
y=
Columns1through8
1.00001.10001.19181.27741.35821.43511.50901.5803
Columns9through11
1.64981.71781.7848
（2）改进Euler算法：
x=
Columns1through7
00.10000.20000.30000.40000.50000.6000
Columns8through11
0.70000.80000.90001.0000
y=
Columns1through7
1.00001.09591.18411.26621.34341.41641.4860
Columns8through11
1.55251.61651.67821.7379
3、Runge－Kutta算法程序；
function[x,y]=Runge_Kutta4（ydot_fun,x0,y0,h,N）
%标准四阶Runge_Kutta公式，其中
%ydot_fun---一阶微分方程的函数
%x0,y0---初始条件
%h---区间步长
%N---区间的个数
%x---Xn构成的向量
%y---Yn构成的向量
x=zeros（1,N+1）;y=zeros（1,N+1）;x
（1）=x0;y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
k1=h*feval（ydot_fun,x（n）,y（n））;
k2=h*feval（ydot_fun,x（n）+1/2*h,y（n）+1/2*k1）;
k3=h*feval（ydot_fun,x（n）+1/2*h,y（n）+1/2*k2）;
k4=h*feval（ydot_fun,x（n）+h,y（n）+k3）;
y（n+1）=y（n）+1/6*（k1+2*k2+2*k3+k4）;
end
4、用Runge－Kutta算法求解P178例题、P285实验任务
（2），计算结果如下（其中表示数值解，表示解析解，结果保留八位有效数字）：
求解P178例题：
x=
00.10000.20000.30000.40000.5000
y=
1.00001.11111.25001.42861.66672.0000
0.05
0.1
0.15
0.2
0.25
-0.0222
-0.0388
-0.0498
-0.0552
-0.0550
-0.0222
-0.0388
-0.0498
-0.0552
-0.0550
0.3
0.35
0.4
0.45
0.5
-0.0493
-0.0380
-0.0213
0.0010
0.0288
-0.0493
-0.0380
-0.0213
0.0010
0.0288
5、P285实验任务
（2）精确解与近似解的图形比较
6、用matlab中的ode45求解P285实验任务
（2）
ans=
Columns1through7
00.00010.00020.00030.00040.00090.0014
0-0.0001-0.0001-0.0002-0.0002-0.0005-0.0007
Columns8through14
0.00190.00240.00490.00740.00990.01250.0250
-0.0010-0.0012-0.0024-0.0037-0.0049-0.0061-0.0118
Columns15through21
0.03750.05000.06250.07500.08750.10000.1125
-0.0172-0.0222-0.0268-0.0312-0.0351-0.0388-0.0420
Columns22through28
0.12500.13750.15000.16250.17500.18750.2000
-0.0450-0.0475-0.0497-0.0516-0.0532-0.0543-0.0552
Columns29through35
0.21250.22500.23750.25000.26250.27500.2875
-0.0556-0.0558-0.0556-0.0550-0.0541-0.0528-0.0512
Columns36through42
0.30000.31250.32500.33750.35000.36250.3750
-0.0493-0.0470-0.0444-0.0414-0.0381-0.0344-0.0304
Columns43through49
0.38750.40000.41250.42500.43750.45000.4625
-0.0260-0.0213-0.0162-0.0108-0.00510.00100.0074
Columns50through53
0.47180.48120.49060.5000
0.01250.01770.02320.0288
五、讨论分析
误差一开始变大，然后变小，最后又变大
六、改进实验建议
可以通过提高保留的位数，使Runge－Kutta算法结果更精确，也可以使数值解和解析解的比较更加精确
实验四常微分方程初值问题数值解法
（这里只提供解答格式请同学自己按照上机实验报告格式填写）
1、饿uler法求解初值问题
function[x,y]=euler_f（ydot_fun,x0,y0,h,N）
%Euler（向前）公式，其中
%ydot_fun---一阶微分方程的函数
%x0,y0---初始条件
%h---区间步长
%N---区间的个数
%x---Xn构成的向量
%y---Yn构成的向量
x=zeros（1,N+1）;y=zeros（1,N+1）;x
（1）=x0;y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
y（n+1）=y（n）+h*feval（ydot_fun,x（n）,y（n））;
end
用欧拉法计算p167例1
>>ydot_fun=inline（'y-2*x./y','x','y'）;
>>[x,y]=euler_f（ydot_fun,0,1,0.1,10）
x=Columns1through11
00.1000000000000000.2000000000000000.3000000000000000.4000000000000000.500000000000000
0.6000000000000000.7000000000000000.8000000000000000.9000000000000001.000000000000000
y=.0000000000000001.1000000000000001.1918181818181821.2774378337147221.3582125995602891.435132918657796
1.5089662535663321.5803382376552171.6497834310477111.7177793478600871.784770832497982
2、改进Euler公式
function[x,y]=euler_r（ydot_fun,x0,y0,h,N）
%改进Euler公式，其中
%ydot
升级会员 