四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材.docx
《四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材.docx（31页珍藏版）》请在冰点文库上搜索。
四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材
实验九欧拉方法
信息与计算科学金融崔振威201002034031
一、实验目的：
1、掌握欧拉算法设计及程序实现
二、实验内容：
1、p364-9.2.4、p386-9.5.6
三、实验要求：
主程序：
欧拉方法（前项）：
function[x,y]=DEEuler（f,a,b,y0,n）;
%f:
一阶常微分方程的一般表达式的右端函数
%a:
自变量的取值下限
%b:
自变量的取值上限
%y0：
函数的初值
%n:
积分的步数
ifnargin<5,n=50;
end
h=（b-a）/n;
x
（1）=a;y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y（i+1）=y（i）+h*feval（f,x（i）,y（i））;
end
formatshort
欧拉方法（后项）：
function[x,y]=BAEuler（f,a,b,y0,n）;
%f:
一阶常微分方程的一般表达式的右端函数
%a:
自变量的取值下限
%b:
自变量的取值上限
%y0：
函数的初值
%n:
积分的步数
ifnargin<5,n=50;
end
h=（b-a）/n;
x
（1）=a;y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y1=y（i）+h*feval（f,x（i）,y（i））;
y（i+1）=y（i）+h*feval（f,x（i+1）,y1）;
end
formatshort
梯形算法：
function[I,step,h2]=CombineTraprl（f,a,b,eps）
%f被积函数
%a,b积分上下限
%eps精度
%I积分结果
%step积分的子区间数
if（nargin==3）
eps=1.0e-4;
end
n=1;
h=（b-a）/2;
I1=0;
I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
whileabs（I2-I1）>eps
n=n+1;
h=（b-a）/n;
I1=I2;
I2=0;
fori=0:
n-1
x=a+h*i;
x1=x+h;
I2=I2+（h/2）*（subs（sym（f）,findsym（sym（f））,x）+subs（sym（f）,findsym（sym（f））,x1））;
end
end
I=I2;
step=n;
h2=（b-a）/n;
改进欧拉方法：
function[x,y]=MoEuler（f,a,b,y0,n）;
%f:
一阶常微分方程的一般表达式的右端函数
%a:
自变量的取值下限
%b:
自变量的取值上限
%y0：
函数的初值
%n:
积分的步数
ifnargin<5,n=50;
end
h=（b-a）/n;
x
（1）=a;y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y1=y（i）+h*feval（f,x（i）,y（i））;
y2=y（i）+h*feval（f,x（i+1）,y1）;
y（i+1）=（y1+y2）/2;
end
formatshort
四阶龙格-库塔法：
functiony=DELGKT4_lungkuta（f,h,a,b,y0,varvec）
%f:
一阶常微分方程的一般表达式的右端函数
%h:
积分步长
%a：
自变量的取值下限
%b:
自变量的取值上限
%varvec：
常微分方程的变量组
formatlong;
N=（b-a）/h;
y=zeros（N+1,1）;
y
（1）=y0;
x=a:
h:
b;
var=findsym（f）;
fori=2:
N+1
K1=Funval（f,varvec,[x（i-1）y（i-1）]）;%Funval为程序所需要的函数
K2=Funval（f,varvec,[x（i-1）+h/2y（i-1）+K1*h/2]）;
K3=Funval（f,varvec,[x（i-1）+h/2y（i-1）+K2*h/2]）;
K4=Funval（f,varvec,[x（i-1）+hy（i-1）+h*K3]）;
y（i）=y（i-1）+h*（K1+2*K2+2*K3+K4）/6;
end
formatshort;
p364-9.2.4
欧拉方法（前项）：
1、
解:
执行20步时：
编写函数文件doty.m如下：
functionf=doty（x,y）;
f=x^2-y;
在Matlab命令窗口输入：
>>[x1,y1]=DEEuler（'doty',0,2,1,20）
可得到结果：
x1=
00.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000
y1=
1.00000.90000.81100.73390.66950.61860.58170.55950.55260.56130.58620.62760.68580.76120.85410.96471.09321.23991.40491.58841.7906
在Matlab命令窗口输入：
>>y3=-exp（-x1）+x1.^2-2*x1+2
求得解析解：
y3=
1.00000.90520.82130.74920.68970.64350.61120.59340.59070.60340.63210.67710.73880.81750.91341.02691.15811.30731.47471.66041.8647
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
执行40步时：
在Matlab命令窗口输入：
>>[x1,y1]=DEEuler（'doty',0,2,1,40）
可得到结果：
x1=
00.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60000.65000.70000.75000.80000.85000.90000.95001.00001.05001.10001.15001.20001.25001.30001.35001.40001.45001.50001.55001.60001.65001.70001.75001.80001.85001.90001.95002.0000
y1=
1.00000.95000.90260.85800.81620.77740.74170.70910.67980.65380.63120.61210.59670.58480.57670.57240.57190.57530.58260.59400.60940.62900.65260.68050.71260.74900.78970.83470.88410.93790.99611.05881.12601.19771.27391.35471.44011.53011.62471.72401.8279
在Matlab命令窗口输入：
>>y3=-exp（-x1）+x1.^2-2*x1+2
求得解析解：
y3=
1.00000.95130.90520.86180.82130.78370.74920.71780.68970.66490.64350.62560.61120.60050.59340.59010.59070.59510.60340.61580.63210.65260.67710.70590.73880.77600.81750.86330.91340.96791.02691.09031.15811.23051.30731.38871.47471.56531.66041.76021.8647
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。
5、
解:
执行20步时：
编写函数文件doty.m如下：
functionf=doty（x,y）;
f=2*x*y^2;
在Matlab命令窗口输入：
>>[x1,y1]=DEEuler（'doty',0,2,1,20）
可得到结果：
x1=
00.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000
y1=
1.0e+176*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00006.4573
在Matlab命令窗口输入：
>>y3=1./（1-x1.^2）
求得解析解：
y3=
1.0e+015*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00004.5036-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
执行40步时：
在Matlab命令窗口输入：
>>[x1,y1]=DEEuler（'doty',0,2,1,40）
可得到结果：
x1=
0.00000.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60000.65000.70000.75000.80000.85000.90000.95001.00001.0500
1.10001.15001.20001.25001.30001.35001.40001.45001.50001.55001.60001.65001.70001.75001.80001.85001.90001.95002.0000
y1=
1.0e+224*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00003.6893InfInfInfInfInfInfInfInfInf在Matlab命令窗口输入：
>>y3=1./（1-x1.^2）
求得解析解：
y3=
1.0e+015*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000-2.2518-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
欧拉方法（后项）：
1、
解:
执行20步时：
编写函数文件doty.m如下：
functionf=doty（x,y）;
f=x^2-y;
在Matlab命令窗口输入：
>>[x1,y1]=BAEuler（'doty',0,2,1,20）
可得到结果：
x1=
0.00000.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000
y1=
1.00000.91100.83290.76650.71270.67190.64490.63230.63450.65200.68520.73450.80030.88290.98251.09951.23411.38641.55671.74521.9520
在Matlab命令窗口输入：
>>y3=-exp（-x1）+x1.^2-2*x1+2
求得解析解：
y3=
1.00000.90520.82130.74920.68970.64350.61120.59340.59070.60340.63210.67710.73880.81750.91341.02691.15811.30731.47471.66041.8647
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
执行40步时：
在Matlab命令窗口输入：
>>[x1,y1]=BAEuler（'doty',0,2,1,40）
可得到结果：
x1=
00.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60000.65000.70000.75000.80000.85000.90000.95001.00001.05001.10001.15001.20001.25001.30001.35001.40001.45001.50001.55001.60001.65001.70001.75001.80001.85001.90001.95002.0000
y1=
1.00000.95260.90790.86580.82670.79040.75720.72720.70030.67680.65660.63990.62680.61720.61140.60920.61090.61640.62580.63920.65660.67800.70350.73320.76710.80520.84750.89420.94511.00051.06021.12431.19291.26601.34351.42561.51221.60341.69921.79961.9046
在Matlab命令窗口输入：
>>y3=-exp（-x1）+x1.^2-2*x1+2
求得解析解：
y3=
1.00000.95130.90520.86180.82130.78370.74920.71780.68970.66490.64350.62560.61120.60050.59340.59010.59070.59510.60340.61580.63210.65260.67710.70590.73880.77600.81750.86330.91340.96791.02691.09031.15811.23051.30731.38871.47471.56531.66041.76021.864
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。
但是后项欧拉要比前项接近解析解。
5、
解:
执行20步时：
编写函数文件doty.m如下：
functionf=doty（x,y）;
f=2*x*y^2;
在Matlab命令窗口输入：
>>[x1,y1]=BAEuler（'doty',0,2,1,20）
可得到结果：
x1=
0.00000.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000
y1=
1.0e+119*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00004.0544InfInfInfInfInfInfInf
在Matlab命令窗口输入：
>>y3=1./（1-x1.^2）
求得解析解：
y3=
1.0e+015*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00004.5036-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
执行40步时：
在Matlab命令窗口输入：
>>[x1,y1]=BAEuler（'doty',0,2,1,40）
可得到结果：
x1=
00.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60000.65000.70000.75000.80000.85000.90000.95001.00001.05001.10001.15001.20001.25001.30001.35001.40001.45001.50001.55001.60001.65001.70001.75001.80001.85001.90001.95002.0000
y1=
1.0e+195*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00001.0834InfInfInfInfInfInfInfInfInfInfInfInfInfInfInfInfInf
在Matlab命令窗口输入：
>>y3=1./（1-x1.^2）
求得解析解：
y3=
1.0e+015*
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000-2.2518-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000
输入：
>>plot（x1,y1,'o'）;holdon
>>plot（x1,y3,'*'）;holdon
可得近似值与解析解的图像比较：
梯形算法：
1、
解：
在matlab窗口中输入：
>>[I,step,h2]=CombineTraprl（'-exp（-x1）+x1^2-2*x1+2',0,2）
可得结果：
I=
1.8032
step=
29
h2=
0.0690
由此可以看出，梯形算法求得的结果比执行20步的欧拉算法要接近解析解，但执行40步的近似解要比梯形算法精确。
5、
解：
在matlab窗口中输入：
>>[I,step,h2]=CombineTraprl（'1/（1-x^2）',0,2）
无法得出结果，程序报
Warning:
Dividebyzero.
错误
改进欧拉方法：
1、
解:
执行20步时：
编写函数文件doty.m如下：
functionf=doty（x,y）;
f=x^2-y;
在Matlab命令窗口输入：
>>[x1,y1]=MoEuler（'doty',0,2,1,20）
可得到结果：
x1=
0.00000.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000
y1=
1.00000.90550.82190.75010.
升级会员 