实验四常微分方程初值问题数值解法.docx
- 文档编号:8699690
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:13
- 大小:79.07KB
实验四常微分方程初值问题数值解法.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_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
用改进欧拉公式计算P167例题1
>>ydot_fun=inline('y-2*x./y','x','y');
>>[x,y]=euler_r(ydot_fun,0,1,0.1,10)
x=Columns1through6
00.1000000000000000.2000000000000000.3000000000000000.4000000000000000.500000000000000
Columns7through11
0.6000000000000000.7000000000000000.8000000000000000.9000000000000001.000000000000000
y=Columns1through6
1.0000000000000001.0959090909090911.184********29971.2662013608757761.3433601514839991.416401928536909
Columns7through11
1.4859556024156691.5525140913261461.6164747827520581.6781663636751861.737867401035414
用欧拉法、改进欧拉法计算P171例题3
ydot_fun=inline('-y+x+1','x','y');
>>[x,y]=euler_f(ydot_fun,0,1,0.1,5)
x=00.1000000000000000.2000000000000000.3000000000000000.4000000000000000.500000000000000
y=1.0000000000000001.0000000000000001.0100000000000001.0290000000000001.0561000000000001.090490000000000
用改进欧拉法计算P171例题3
ydot_fun=inline('-y+x+1','x','y');
>>[x,y]=euler_r(ydot_fun,0,1,0.1,5)
x=00.1000000000000000.2000000000000000.3000000000000000.4000000000000000.500000000000000
y=1.0000000000000001.0050000000000001.0190250000000001.0412176250000001.0708019506250001.107075765315625
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
用四阶龙格库塔计算P178例题
>>ydot_fun=inline('y^2','x','y');
>>[x,y]=Runge_Kutta4(ydot_fun,0,1,0.1,5)
x=00.1000000000000000.2000000000000000.3000000000000000.4000000000000000.500000000000000
y=1.0000000000000001.1111104900521941.2499979920470151.4285661863014451.6666532572503231.999963258950669
用龙格库塔方法计算P285实验任务
(2)
>>ydot_fun=inline('(-y+x^2+4*x-1)./2','x','y');
>>[x,y]=Runge_Kutta4(ydot_fun,0,0,0.05,10)
x=Columns1through6
00.0500000000000000.1000000000000000.150********00000.2000000000000000.250000000000000
Columns7through11
0.3000000000000000.3500000000000000.4000000000000000.4500000000000000.500000000000000
y=Columns1through6
.022*********
Columns7through11
-0.049292018554665-0.038042973450910-0.0212692404030080.0010162259975860.028800791009418
>>xx=[0:
0.05:
0.5];
>>z=exp(-xx./2)+xx.^2-1;
>>holdon
>>plot(x,y,'o');plot(xx,z,'*')
>>holdoff
其图形如图一所示
图一精确解与龙格库塔计算结果比较
5、用Matlab提供的ode45求解P285实验任务
(2)
>>ydot_fun=inline('(-y+x^2+4*x-1)./2','x','y');
>>[x,y]=ode45(ydot_fun,[0,0.5],0);
>>[x';y']
>>plot(x,y,’*’)
ans=Columns1through6
00.0001004754572600.0002009509145210.0003014263717810.0004019018290420.0009042791153430-0.000050226371419-0.000100430028501-0.000150610971371-0.000200769200158-0.000451219637267
0.0014066564016450.0019090336879470.0024114109742490.0049232974057590.0074351838372680.009947070268778
-0.000701102241287-0.000950417028058-0.001199164013417-0.002434382472993-0.003655408270323-0.004862243380400
0.0124589567002880.0249589567002880.0374589567002880.0499589567002880.0624589567002880.074958956700288
-0.006054889775763-0.011778983079828-0.017151998201403-0.022174175468426-0.026845753781789-0.031166970574532
.0874********
.0351********
0.1624589567002880.1749589567002880.187********02880.199********02880.2124589567002880.224958956700288
.0516********
0.2374589567002880.2499589567002880.2624589567002880.2749589567002880.2874589567002880.299958956700287
.0555********
0.3124589567002880.3249589567002880.3374589567002870.3499589567002870.3624589567002870.374958956700288
.0470********
0.3874589567002880.3999589567002870.4124589567002870.4249589567002870.4374589567002870.449958956700287
.025*********
0.4624589567002870.4718442175252160.4812294783501440.4906147391750720.500000000000000
0.0074256275471370.0124791589333560.0177262495359700.0231668179352700.028800783076419
图二用ode45求解得到的近似解散点图
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验四 常微分方程初值问题数值解法 实验 微分方程 初值问题 数值 解法