实验四常微分方程初值问题数值解法.docx

实验四常微分方程初值问题数值解法.docx

《实验四常微分方程初值问题数值解法.docx》由会员分享，可在线阅读，更多相关《实验四常微分方程初值问题数值解法.docx（17页珍藏版）》请在冰豆网上搜索。

实验四常微分方程初值问题数值解法

科学技术学院

实验报告

课程名称数值分析

实验项目常微分方程初值问题数值解法

专业班级12.数学与应用数学（师）叶楚欣学号2012214103

指导教师黄国顺成绩日期

一.实验目的

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算法程序；

2、用Euler算法程序、改进Euler算法求解P167例题1的运行结果；

3、Runge－Kutta算法程序；

4、用Runge－Kutta算法求解P178例题、P285实验任务

（2），计算结果如下（其中

表示数值解，

表示解析解，结果保留八位有效数字）：

0.05

0.1

0.15

0.2

0.25

-0.02219009

-0.03877057

-0.04975651

-0.05516258

-0.05500309

-0.02219009

-0.03877058

-0.04975651

-0.05516258

-0.05500310

0.3

0.35

0.4

0.45

0.5

-0.04929202

-0.03804297

-0.02126924

0.00101623

0.02880079

-0.04929202

-0.03804298

-0.02126925

0.00101622

0.02880078

5、P285实验任务

（2）精确解与近似解的图形比较

（图片贴到此处）

6、用matlab中的ode45求解P285实验任务

（2）

MATLABdesktopkeyboardshortcuts,suchasCtrl+S,arenowcustomizable.

Inaddition,manykeyboardshortcutshavechangedforimprovedconsistency

acrossthedesktop.

Tocustomizekeyboardshortcuts,usePreferences.Fromthere,youcanalso

restorepreviousdefaultsettingsbyselecting"R2009aWindowsDefaultSet"

fromtheactivesettingsdrop-downlist.Formoreinformation,seeHelp.

Clickhereifyoudonotwanttoseethismessageagain.

>>ydot_fun=inline（'y-2*x./y','x','y'）;

>>[x,y]=euler_f（ydot_fun,0,1,0.1,10）

x=

Columns1through4

00.00000.00000.0000

Columns5through8

0.00000.00000.00000.0000

Columns9through11

0.00000.00001.0000

y=

Columns1through4

1.00001.00001.81821.4722

Columns5through8

1.02891.77961.63321.5217

Columns9through11

1.77111.00871.7982

>>ydot_fun=inline（'y-2*x./y','x','y'）;

>>[x,y]=euler_r（ydot_fun,0,1,0.1,10）

x=

Columns1through4

00.00000.00000.0000

Columns5through8

0.00000.00000.00000.0000

Columns9through11

0.00000.00001.0000

y=

Columns1through4

1.00001.90911.29971.5776

Columns5through8

1.39991.69091.56691.6146

Columns9through11

1.20581.51861.5414

>>ydot_fun=inline（'-y+x+1','x','y'）;

>>[x,y]=euler_f（ydot_fun,0,1,0.1,5）

x=

Columns1through4

00.00000.00000.0000

Columns5through6

0.00000.0000

y=

Columns1through4

1.00001.00001.00001.0000

Columns5through6

1.00001.0000

>>ydot_fun=inline（'-y+x+1','x','y'）;

>>[x,y]=euler_r（ydot_fun,0,1,0.1,5）

x=

Columns1through4

00.00000.00000.0000

Columns5through6

0.00000.0000

y=

Columns1through4

1.00001.00001.00001.0000

Columns5through6

1.50001.5625

>>ydot_fun=inline（'y^2','x','y'）;

>>[x,y]=Runge_Kutta4（ydot_fun,0,1,0.1,5）

x=

Columns1through4

00.00000.00000.0000

Columns5through6

0.00000.0000

y=

Columns1through4

1.00001.21941.70151.1445

Columns5through6

1.03231.0669

>>ydot_fun=inline（'（-y+x^2+4*x-1）./2','x','y'）;

>>[x,y]=Runge_Kutta4（ydot_fun,0,0,0.05,10）

x=

Columns1through4

00.00000.00000.0000

Columns5through8

0.00000.00000.00000.0000

Columns9through11

0.00000.00000.0000

y=

Columns1through4

0-0.6823-0.3692-0.8554

Columns5through8

-0.6671-0.5772-0.4665-0.0910

Columns9through11

-0.30080.75860.9418

ydot_fun=inline（'（-y+x^2+4*x-1）./2','x','y'）;

[x,y]=Runge_Kutta4（ydot_fun,0,0,0.05,10）;sprintf（'%.8f,%.8f\n',x,y）

ans=

0.00000000,0.05000000

0.10000000,0.15000000

0.20000000,0.25000000

0.30000000,0.35000000

0.40000000,0.45000000

0.50000000,0.00000000

-0.02219009,-0.03877057

-0.04975651,-0.05516258

-0.05500309,-0.04929202

-0.03804297,-0.02126924

0.00101623,0.02880079

xx=[0:

0.05:

0.5];

z=exp（-xx./2）+xx.^2-1;sprintf（'%.8f,%.8f\n',xx,z）

ans=

0.00000000,0.05000000

0.10000000,0.15000000

0.20000000,0.25000000

0.30000000,0.35000000

0.40000000,0.45000000

0.50000000,0.00000000

-0.02219009,-0.03877058

-0.04975651,-0.05516258

-0.05500310,-0.04929202

-0.03804298,-0.02126925

0.00101622,0.02880078

>>xx=[0:

0.05:

0.5];

z=exp（-xx./2）+xx.^2-1;

holdon

>>plot（x,y,'o'）;plot（xx,z,'*'）;

>>holdoff

>>ydot_fun=inline（'（-y+x^2+4*x-1）./2','x','y'）;

>>[x,y]=ode45（ydot_fun,[0,0.5],0）;

>>[x';y']

ans=

Columns1through4

00.72600.45210.1781

0-0.1419-0.8501-0.1371

Columns5through8

0.90420.53430.16450.7947

-0.0158-0.7267-0.1287-0.8058

Columns9through12

0.42490.57590.72680.8778

-0.3417-0.2993-0.0323-0.0400

Columns13through16

0.02880.02880.02880.0288

-0.5763-0.9828-0.1403-0.8426

Columns17through20

0.02880.02880.02880.0288

-0.1789-0.4532-0.3373-0.1758

Columns21through24

0.02880.02880.02880.0288

-0.2876-0.8809-0.2155-0.8392

Columns25through28

0.02880.02880.02880.0288

-0.8617-0.8665-0.0809-0.5467

Columns29through32

0.02880.02880.02880.0288

-0.3663-0.6119-0.5031-0.5134

Columns33through36

0.02880.02880.02880.0287

-0.5891-0.0564-0.8055-0.3352

Columns37through40

0.02880.02880.02870.0287

-0.9452-0.6404-0.3229-0.7744

Columns41through44

0.02870.02880.02880.0287

-0.8238-0.2486-0.9740-0.9953

Columns45through48

0.02870.02870.02870.0287

-0.5211-0.8719-0.68590.2165

Columns49through52

0.02870.52160.01440.5072

0.71370.33560.59700.5270

Column53

0.0000

0.6419

>>plot（x,y,'*'）

>>

五、讨论分析

（通过最后绘制的图形直观分析近似解与准确解之间的误差，自己补充）

从图中可以看到4阶的Runge-Kutta画出来的点基本和原函数的对应点相重合，精度之高，从数据来看，起码有9位有效数字。

六、改进实验建议

（自己补充）

可以通过输出数据的精度控制函数

实验四常微分方程初值问题数值解法

（这里只提供解答格式请同学自己按照上机实验报告格式填写）

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.00000.00000.00000.00000.0000

0.00000.00000.00000.00001.0000

y=.00001.00001.81821.47221.02891.7796

1.63321.52171.77111.00871.7982

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.00000.00000.00000.00000.0000

Columns7through11

0.00000.00000.00000.00001.0000

y=Columns1through6

1.00001.90911.29971.57761.39991.6909

Columns7through11

1.56691.61461.20581.51861.5414

用欧拉法、改进欧拉法计算P171例题3

ydot_fun=inline（'-y+x+1','x','y'）;

>>[x,y]=euler_f（ydot_fun,0,1,0.1,5）

x=00.00000.00000.00000.00000.0000

y=1.00001.00001.00001.00001.00001.0000

用改进欧拉法计算P171例题3

ydot_fun=inline（'-y+x+1','x','y'）;

>>[x,y]=euler_r（ydot_fun,0,1,0.1,5）

x=00.00000.00000.00000.00000.0000

y=1.00001.00001.00001.00001.50001.5625

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.00000.00000.00000.00000.0000

y=1.00001.21941.70151.14451.03231.0669

用龙格库塔方法计算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.00000.00000.00000.00000.0000

Columns7through11

0.00000.00000.00000.00000.0000

y=Columns1through6

0-0.6823-0.3692-0.8554-0.6671-0.5772

Columns7through11

-0.4665-0.0910-0.30080.75860.9418

>>xx=[0:

0.05:

0.5];

>>z=exp（-xx./2）+xx.^2-1;

>>holdon

>>plot（x,y,'o'）;plot（xx,z,'*'）

>>holdoff

其图形如图一所示

图一精确解与龙格库塔计算结果比较

4、用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.72600.45210.17810.90420.53430-0.1419-0.8501-0.1371-0.0158-0.7267

0.16450.79470.42490.57590.72680.8778

-0.1287-0.8058-0.3417-0.2993-0.0323-0.0400

0.02880.02880.02880.02880.02880.0288

-0.5763-0.9828-0.1403-0.8426-0.1789-0.4532

0.02880.02880.02880.02880.02880.0288

-0.3373-0.1758-0.2876-0.8809-0.2155-0.8392

0.02880.02880.02880.02880.02880.0288

-0.8617-0.8665-0.0809-0.5467-0.3663-0.6119

0.02880.02880.02880.02880.02880.0287

-0.5031-0.5134-0.5891-0.0564-0.8055-0.3352

0.02880.02880.02870.02870.02870.0288

-0.9452-0.6404-0.3229-0.7744-0.8238-0.2486

0.02880.02870.02870.02870.02870.0287

-0.9740-0.9953-0.5211-0.8719-0.68590.2165

0.02870.52160.01440.50720.0000

0.71370.33560.59700.52700.6419

图二用ode45求解得到的近似解散点图