龙格库塔求解微分方程Word格式.docx
- 文档编号:16510183
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:17
- 大小:477.04KB
龙格库塔求解微分方程Word格式.docx
《龙格库塔求解微分方程Word格式.docx》由会员分享,可在线阅读,更多相关《龙格库塔求解微分方程Word格式.docx(17页珍藏版)》请在冰豆网上搜索。
2经典龙格-库塔算法
公式(15张)
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
由于此
算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
对于一阶精度的欧拉公式有:
yi+1=yi+h*K1
K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均
斜率K*的近似值,那么就会得到二阶精度的改进欧拉公式:
yi+1=yi+h*(K1+K2)/2
K2=f(xi+h,yi+h*K1)
依次类推,如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km并
用他们的加权平均数作为平均斜率K*的近似值,显然能构造出具有很高精度的高
阶计算公式。
经数学推导、求解,可以得出四阶龙格-库塔公式,也就是在工程中应用广泛的经典龙格-库塔算法:
yi+1=yi+h*(K1+2*K2+2*K3+K4)/6
K2=f(xi+h/2,yi+h*K1/2)
K3=f(xi+h/2,yi+h*K2/2)
K4=f(xi+h,yi+h*K3)
3MATLAB算法的实现
3.1程序代码:
function[tpyp]=RK45_n_h(dydt,tspan,yO,h)formatIong;
tf=tspan
(2);
n=length(t);
%
whilet(n)<
tf
t(n+1)=t(n)+hh;
n=n+1;
ift(n)>
t(n)=tf;
break;
t=(tspan
(1):
h:
tspan
(2))'
;
%B
因为在t=(tspan
(1):
tspan
(2))'
累加过程中,
可能出现大数吃小数的可能,这样可以确保积到上限
An依然指向积分上限
Bt依然等间隔的存储着积分上下限
end
tt=tspan
(1);
%yp(np,:
)用来记录函数值随时间的变化,%为了作图的精确方便,就以步长为间隔.
y(1,:
)=y0;
np=1;
tp(np)=tt;
yp(np,:
)=y(1,:
);
i=1;
while
(1)
tend=t(np+1);
hh=t(np+1)-t(np);
ifhh>
h,hh=h;
endwhile
(1)
iftt+hh>
tend,hh=tend-tt;
end
k1=dydt(tt,y(i,:
))'
k2=dydt(tt+hh/2,y(i,:
)+k1.*hh./2)'
k3=dydt(tt+hh/2,y(i,:
)+k2.*hh./2)'
k4=dydt(tt+hh,y(i,:
)+k3*hh)'
y(i+1,:
)=y(i,:
)+(k1+2*(k2+k3)+k4)/6*hh;
i=i+1;
tt=t(i);
iftt>
=tend,break,end
np=np+1;
tp(np)=tt;
yp(np,:
%溢出时跳出循环
ifmax(y(np,:
))>
1e200
=tf,break,end
formatshort
3.2应用力证
function[tpyp]=RK45_n_h(dydt,tspan,yO,h)
dydt是引用方程组,tspan是积分上下限;
y0是初值;
h为步长.
函数调用格式:
[AB]=RK45_n_h(@equition,[0,10],[0.1,0.1],0.01)
函数返回结果是两个向量tp和yp,tp是从0开始,以0.01为步长倒10的向量;
yp是对应的函数值向量,有两个分量.初始值[0.10.1].
输入命令:
plot(B(:
1),B(:
2)),grid;
得到实时解的平面相轨迹
从图中反映出其解具有周期性,轨迹构成稳定的极限环.
如果初值远离圆环X2+y2=1,如[10,10],解又是怎样的?
[AB]=RK45_n_h(@equition,[0,10],[10,10],0.01);
plot(B(:
2));
grid
解仍然是稳定的极限环
为验证其准确性,可取时间跨度大一些[0,50],重复上述操作,事实上还可以更直观的观察解的变化:
plot(A,B(:
1),A,B(:
2)),grid
可以看出随时间推移,x
(1),x
(2)两个分量都是简谐振动,相位相差
小结
以上就是基于MATLAB平台数值方法解常微分方程的大致过程
其他二维自治微分方程组的解法举例
4・1・1・a
dx
dt
dy
X(0)
2k-xy
-y+xy
=xOfy{0)=yO
分别取不同的初值[1,0.2];
[1,0.5];
[10.8],[1,1.1];
[1,1.2]
axis([O606]);
xlabel('
x轴'
ylabel('
y轴'
[t,x]=ode23(ode913,[06],[10.2],option);
[t,x]=ode23(ode913,[06],[10.5],option);
[t,x]=ode23(ode913,[06],[10.8],option);
[t,x]=ode23(ode913,[06],[11.1],option);
[t,x]=ode23(ode913,[06],[11.2],option);
holdon
legend('
初值为[1,0.2]'
'
初值为[1,0.5]'
初值为[1,0.8]'
初值为[1,1.1]'
初值为[1,1.2]'
text(4,0.5,'
text(3.5,2,'
text(2.2,2.5,'
text(1.3,3,'
text(0.9,2,'
得到:
S,
4
初佰为[■!
(7]值为11.0.51
为仲]
5
百一初值为卩,n.2]
初情対[1』5]
A初11为卩,0.3]
『初11^)11,1.1]
A初侑为ru.2]
5U
J
SloeI
廿由
不同初值,有不同的极限环系统是稳定的
4・1・1・b求如下系统的极限环,并判断其稳定性
■-I,2.2\1/2/2,2八2,
IX=x(x+y)(x+y-1)+y
(L/2,2\1/2/2,2八2
=y(x+y)(X+y-1)-x
输入:
[AB]=RK45_n_h(@BWD,[0,5],[0.001,0.00001],0.001);
plot(A,B(:
BFigure1
FikEditVia神InsertTaoleDefktapWindowHelp
DWHSl釘愆0霍®
宦□SIJ□
■/
D.OB
0.07
D.05
UlUT
2J
[AB]=RK45_n_h(@BWD,[0,10],[0.01,0.1],0.001);
plot(A,B(:
K1FigureI
Fil专EditVi日旳InsertTooItDesktopWindowHelp
Qge毎丨仿迤直■o®
运丨tJE■□
可见初值不同,只决定了开始是的轨迹,但在tT+=^的时候函数值也是趋于无穷的,(尽管在从R=1的圆内趋于圆的速度远小于从圆外趋于无穷的速度)是不稳定的
i.范德波尔方程(解高阶常微分方程)
X+E&
2T)X+x=0
将范德波尔方程化为状态方程
阡y
rI2
$Y(1—X)y—X
取7时:
编写M.函数
%定义输入,输出变量和函数文件名
functiondy=verderpoI(t,y)dy-[y
(2);
7*(1-y
(1)A2)*y
(2)-y
(1)]
ISI回23
QFigure1****
FileEditVinwIreertTddIbDesktopWindow/HalpD®
H爭型理I^HO
15
+特接函数
*一阶异函数
10
4.1.3吸引子
吸引子能量耗散系统最终收缩到的一种定常状态。
这是一个动力系统在t时所呈现的与时
间无关的定态,并且不管选取什么样的初始值其终值的定态只有一个,也就是说终值与初始值无关。
这类吸引子也称平庸吸引子。
如:
阻尼单摆有不动点吸引子,范德玻耳方程有极限环吸引子,等等。
奇怪吸引子相对于平庸吸引子而言,它们的特点之一是终态值与初始值密切相关,或者说对初始值具有极端敏感性;
初始取值的细微差别可能会导致完全不同的结果,这时的吸引子毫无周期可言,即所谓混沌。
4・1・3・aLorenz模型状态方程
rx
(X-y)
xz
Idz
-bz+xy
在洛伦兹方程中,取参数s=10,b=8/3,随参数r增加,出现一次新分岔-霍夫分岔取r=28时计算的结果如下。
dxl/dt=-8/3*xl+x2*x3dx2/dt=-10*x2+10*x3dx3/dt=-xl*x2+28*x2-x3
创建M.函数.
functionxdot=lorenzq(t,x)xdot=[-8/3*x
(1)+x
(2)*x(3);
-10*x
(2)+10*x(3);
-x
(1)*x
(2)+28*x
(2)-x(3)];
在命令窗口输入
[AB]=RK45_n_h(@lorenzq,[0,50],[0,0,1e-10],0.001);
plot(A,B),grid
将X三个分量反映到的三维空间中得到下图
plot3(B(:
2),B(:
3));
不同时刻观察
50
t=5
t=50
30
20
to
JO
-20
ID
从不同角度观察
>
t=60
-5
1
Art
L
]
丄
zf
■
E嚟巍熒
牛x-y平面
奇怪吸引子的最重要特征是对初值的敏感性,初始相互靠近的两条轨线将按指数式规律分离。
但在有限空间中如何保持这样的指数式分离状态?
洛伦兹吸引子有两个不稳定平衡点,因此复杂的相轨线可以随机地在两个中心之间行走。
是否只有一个平衡点的奇怪吸引子呢?
4.1.3.bRdssler方程组
取参数a=b=0.2,c=5.7时计算得罗斯勒吸引子图象
functiondydt=Rossler(t,x)dydt=[-(x
(2)+x(3));
x
(1)+0.2*x
(2);
0.2+(x
(1)-5.7)*x(3)]
[AB]=RK45_n_h(@Rossler,[0,500],[0,0,1e-10],0.01);
■aFpfcjre1
11=1;
LsJ
E3
1File-EcJi-tViewInje-rtToolsOe^k-top>
#Vindc5wHelp
D]1电包心⑥幻u1^pa
3)):
只有一个平衡点的奇怪吸引子
'
r-
J■
T_-
._JSdBa」qa・\Eia-~・・i
■■I■
y宀u、y7:
:
3
■「金0:
■='
b-F--——.
I
参数:
a=b=0.2。
当c=2.6时,相轨线是简单单周期的极限环,其功率谱为系统的基频f(~16Hz)及其谐波;
当c=3.5时,得二周期运动相轨线,其功率谱为f/2、f及其谐波;
当c=4.1时,为四周期的极限环,功率谱为f/4、f/2、f及其谐波;
当c=4.18时,为八周期的极限环,功率谱为f/&
f/4、f/2、f及谐波。
可见随c增加,存在一系列时倍周期分岔,直到倍周期积累点。
过积累点后,随C增加,轨线展宽开来,相近相轨线合并形成宽阔相轨道。
在c=4.21时,相轨线仍是八周期极限环,在功率谱上除f/&
f/4、f/2、f及其谐波
外,还有f/16的弱峰。
在c=4.6时,相轨线成了一条粗大的环线,表明运动已无任何周期,功率谱是在很
高的噪声背景谱上存在着频率f及其谐波的尖峰。
后者奇怪吸引子功率谱的特征。
总结
1优点:
(1)经典龙格-库塔精度较高,在数值计算中广泛应用,而基于MATLAB图形可视化编程使其在分析二维自制系统稳定性问题中,应用方便.结果直观.
(2)该程序的可靠性很好,计算结果稳定,可读性好,在某些算法上做了相应优化,,占用内存较小,使得程序在执行速度上有了一定提高
(最大占用内存:
55M用时最长时:
6分钟.与具体问题(方程)积分区间步长有
关).
(3)精度可调,设定不同的步长,可得到不同精度的结果,用时也不一样.
2不足:
(1)不能实现变步长的有效算法.要求有很高精度时,该算法需人为尝试不同的步长,浪
费时间.
(2)解一阶方程时,该算法存在严重浪费内存的问题.建议用下面的算法(可实现变步长的算法)
functionodeRK45_h1(a,b,yO,n,eps)y1=y0;
h=(b-a)/n;
x=a:
b;
y=y0*ones(1,n+1);
forj=2:
n+1
k1=f(x(j-1),y(j-1));
k2=f(x(j-1)+h/2,y(j-1)+h/2*k1);
k3=f(x(j-1)+h/2,y(j-1)+h/2*k2);
k4=f(x(j-1)+h,y(j-1)+h*k3);
y(j)=y(j-1)+h/6*(k1+k4)+h/3*(k2+k3);
%fork=1:
n+1;
%fprintf('
x[%d]=%.7f\ty[%d]=%.7f\n'
k-1,x(k),k-1,y(k));
%end
whileabs(y(j)-y1)>
15*eps
y1=y(j);
n=n*2;
fork=1:
fprintf('
%有待于将两种方法结合起来使用
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 龙格库塔 求解 微分方程