常微分方程初值问题数值解法的比较.docx
- 文档编号:828871
- 上传时间:2022-10-13
- 格式:DOCX
- 页数:20
- 大小:222.69KB
常微分方程初值问题数值解法的比较.docx
《常微分方程初值问题数值解法的比较.docx》由会员分享,可在线阅读,更多相关《常微分方程初值问题数值解法的比较.docx(20页珍藏版)》请在冰豆网上搜索。
常微分方程初值问题数值解法的比较
常微分方程初值问题数值解法的比较
数值计算实践—课程设计报告
课题名称
常微分方程初值问题数值解法的比较
完成时间
2013-1-17
姓名
班级
学号
成绩
一.实验目的及内容
1实验目的:
(1)了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究;
(2)熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析;
(3)比较欧拉法、改进欧拉法及龙格-库塔法,能够选择合适的方法进行问题的研究计算;
2实验内容:
求微分方程(欧拉法求解)
求微分方程(改进欧拉法求解)
求微分方程(龙格-库塔求解)
根据实验的结果进行分析,了解一般方法的的优缺点,稳定性,收敛性以及截断误差的分析,针对相应问题拿出有效方法得出最优的结果。
二.相关背景知识介绍以及初值问题稳定性的研究:
在科学与工程问题中,常微分方程表述物理量的变化规律,应用非常广泛,比如,天体运动的轨迹,机器人控制,化学反应过程的描述和控制以及电路瞬态过程分析等。
这些问题中要求解随时间变化的物理量,即位置函数表示时间,而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。
考虑一阶常微分方程的初值问题,如果存在实数使得则称f关于y满足利普希茨条件,L称为利普希茨常数。
对于常微分方程初值问题,考虑初值的扰动是问题的解发生偏差的情形。
若时的偏差被控制在有界范围内,则称该初值问题是稳定的,否则该初值问题不稳定的。
特别地,若时的偏差收敛于零,则称该初值问题是渐进稳定的。
对于初值问题稳定性的研究,易知其准确解为,假定初值经过扰动后变为,对于扰动后的解为因此带来的扰动误差为,因此考虑时的值,它取决于。
易知,若,则原问题是稳定的;若,则原问题是不稳定的;若,则原问题是渐进稳定的。
实际遇到的大多数常微分方程初值问题都是稳定的,因此在后面的讨论数值解法时这常常是默认条件。
1.欧拉法:
依据:
积分曲线上一点的切线斜率等于函数值。
方法:
推进法,初始点出发,依照方向场在改点的方向推进到
向前欧拉法的得到:
(1)将在处泰勒展开取h的线性部分,得
(2)将初值问题中得导数用向前差商来代替有,因此
(3)将两边同时对x的区间上积分对右端用左矩形公式得,此方法称向前欧拉法,也叫显示欧拉法。
(4)对右端用右矩形公式得,也叫隐式欧拉法。
误差分析:
1.称为计算时的局部截断误差;
2.如果数值方法的局部截断误差为,那么称这种数值方法的阶数是p,其实p为非负整数。
通常情况下,步长h越小,p越高,则局部截断误差越小,计算精
初泰勒展开有
则有可见欧拉方法是一阶方法,精度不是很高。
2.改进欧拉方法:
梯形公式:
对右端用梯形公式得+显然梯形公式是隐式公式。
改进欧拉公式:
先用欧拉公式求的一个初步的近似值,成为预测值,预测值的精度可能达不到要求,在用梯形公式将他校正一次,记为,这个结果成为校正值。
预测:
校正:
误差分析:
记为改进欧拉公式在处的截断误差,
记因此有
,表示在出的局部截断误差。
由此得,梯形公式的局部截断误差为,因此改进欧拉的截断误差为,可见改进欧拉的方法是二阶方法,改进欧拉方法优于欧拉方法。
3.龙格—库塔法:
根据拉格朗日微分中值定理,,记得到,这样,只要给出一种计算的算法,就可以得到相应的计算公式。
欧拉公式可以写为
改进欧拉公式可以写成因此推出一般的推出广式
称为p阶龙格-库塔方法,简称p阶R-K方法。
因为这里的均为常数。
因为给定的系数不唯一,因此这里的常数有无穷多个解,下面是特殊情况下和一般情况下得结果
(一阶龙格-库塔)当r=1时,这就是欧拉法。
(二阶龙格-库塔)当r=2时,,这就是改进欧拉法。
三阶和四阶龙格-库塔也只是在一般情况下得结果。
三阶
四阶
其局部截断误差是:
三.程序代码
欧拉法代码如下:
(1)向前欧拉法:
functiony=Euler1(fun,x0,y0,xN,N)
%fun为一阶微分方程,x0,y0为初始条件,xN为取值范围的一个端点,N为区间个数
x=zeros(1,N+1);
x=zeros(1,N+1);
x
(1)=x0;
y
(1)=y0;
h=(xN-x0)/N;%h为区间步长
for(n=1:
N)
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));%根据向前欧拉公式计算y值
end
T=[x',y']
(2)向后欧拉法:
functiony=Euler2(fun,x0,y0,xN,N)
%fun为一阶微分方程,x0,y0为初始条件,xN为取值范围的一个端点,N为区间个数x=zeros(1,N+1);
x=zeros(1,N+1);
x
(1)=x0;
y
(1)=y0;
h=(xN-x0)/N;%h为区间步长
for(n=1:
N)
x(n+1)=x(n)+h;
z0=y(n)+h*feval(fun,x(n),y(n));
for(k=1:
3)
z1=y(n)+h*feval(fun,x(n+1),z0);
if(abs(z1-z0)<1e-3)
break;
end
z0=z1;
end
y(n+1)=z1;%%根据向后欧拉公式计算y值
end
T=[x',y']
改进欧拉代码如下:
function[x,y]=Gaijineuler(f,x0,y0,xZ,h)
%f为一阶微分方程的函数;x0,y0为初始条件;xZ为取值范围的一个端点,h为区间步长
n=fix((xZ-x0)/h);%计算分点数
y
(1)=y0;
x
(1)=x0;
fori=1:
n
x(i+1)=x0+i*h;
yp=y(i)+h*feval(f,x(i),y(i));
yc=y(i)+h*feval(f,x(i+1),yp);
y(i+1)=(yp+yc)/2;%根据改进欧拉公式计算结果
end
x=x';
y=y';
1.
3.龙格-库塔代码如下:
(三阶龙格-库塔)
functionR=Longgekuta3(f,a,b,aZ,h)
%a,b为端点,h为步长,aZ为初值
n=(b-a)/h;
T=zeros(1,n+1);%定义向量
Y=zeros(1,n+1);
T=a:
h:
b;%计算各个分点
Y
(1)=aZ;%初值赋予
forj=1:
n
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);
k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h);
Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6;%根据公式计算
end
R=[T'Y'];
(四阶龙格-库塔)
functionR=Longgekuta4(f,a,b,aZ,h)
%a,b为端点,h为步长,aZ为初值
n=(b-a)/h;
T=zeros(1,n+1);%定义向量
Y=zeros(1,n+1);
T=a:
h:
b;%计算各个分点
Y
(1)=aZ;%初值赋予
forj=1:
n
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);
k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);
k4=feval(f,T(j)+h,Y(j)+k3*h);
Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6;%根据公式计算
end
R=[T'Y'];
四.数值结果:
输入:
定义M文件并保存ffx.m
Euler1('ffx',0,1,1,10)
结果:
T=
01.0000
0.10001.1000
0.20001.2100
0.30001.3310
0.40001.4641
0.50001.6105
0.60001.7716
0.70001.9487
0.80002.1436
0.90002.3579
1.00002.5937
ans=
Columns1through7
1.00001.10001.21001.33101.46411.61051.7716
Columns8through11
1.94872.14362.35792.5937
向前欧拉公式结果:
输入:
Euler2('ffx',0,1,1,10)
结果T=
01.0000
0.10001.1110
0.20001.2344
0.30001.3716
0.40001.5240
0.50001.6933
0.60001.8814
0.70002.0904
0.80002.3227
0.90002.5807
1.00002.8674
ans=
Columns1through7
1.00001.11101.23441.37161.52401.69331.8814
Columns8through11
2.09042.32272.58072.8674
改进欧拉公式结果:
输入:
[x,y]=Gaijineuler('f',0,1,1,0.1)
结果:
x=
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
y=
1.0000
1.0959
1.1841
1.2662
1.3434
1.4164
1.4860
1.5525
1.6165
1.6782
1.7379
龙格-库塔计算结果:
(三阶)输入:
Longgekuta3('ff',0,1,1,0.1)
结果:
ans=
Longgekuta3('ff',0,1,1,0.1)
ans=
01.0000
0.10001.1048
0.20001.2188
0.30001.3411
0.40001.4711
0.50001.6082
0.60001.7520
0.70001.9021
0.80002.0580
0.90002.2195
1.00002.3863
(四阶)输入:
Longgekuta4('ff',0,1,1,0.1)
结果:
ans=
01.0000
0.10001.1048
0.20001.2188
0.30001.3411
0.40001.4711
0.50001.6082
0.60001.7520
0.70001.9021
0.80002.0580
0.90002.2195
1.00002.3863
五.计算结果分析:
方法
显示欧拉
简单
精度低
隐式欧拉
稳定性最好
精度低,计算量大
梯形公式
精度提高
计算量大
中点公式
精度提高,显式
多一个初值,可能影响精度
1.收敛性:
若某算法对于任意固定的x=xi=x0+ih,当h0(同时i)时有yiy(xi),则称该算法是收敛的。
以下讨论的都是单步法(指在计算时只用到它前一步的
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 初值问题 数值 解法 比较