matlab在数值分析中的应用runge-kutta

《matlab在数值分析中的应用runge-kutta》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、第七章微分方程问题的解法微分方程的解析解方法微分方程问题的数值解法微分方程问题算法概述四阶定步长Runge-Kutta算法及MATLAB实现一阶微分方程组的数值解微分方程转换特殊微分方程的数值解7.1微分方程的解析解方法格式：y=dsolve(f1,f2,…,fm)格式：指明自变量y=dsolve(f1,f2,…,fm,’x’)fi即可以描述微分方程，又可描述初始条件或边界条件。如：描述微分方程时描述条件时例：>>symst;u=exp(-5*t)*cos(2*t+1)+5;>>uu=5*diff

2、(u,t,2)+4*diff(u,t)+2*uuu=87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10>>symsty;>>y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])>>y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...'87*exp(-5*t)*cos(

3、2*t+1)+92*exp(-5*t)*sin(2*t+1)...+10'],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')分别处理系数，如：>>[n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))]ans=-8704185%rat()最接近有理数的分数判断误差：>>vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)ans=.1147319758647

4、90922564144636e-4>>y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+...10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5');如果用推导的方法求Ci的值，每个系数的解析解至少要写出10数行,故可采用有理式近似的方式表示.>>vpa(y,10)%有理近似值ans=1.19636183

5、9*exp(-5.*t)+.4166666667-.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e-1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(-5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(-3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(-1.*t)+442590.9059*exp(-4.

6、*t)例：>>symstx>>x=dsolve('Dx=x*(1-x^2)')x=[1/(1+exp(-2*t)*C1)^(1/2)][-1/(1+exp(-2*t)*C1)^(1/2)]>>symstx;x=dsolve('Dx=x*(1-x^2)+1')Warning:Explicitsolutioncouldnotbefound;implicitsolutionreturned.>InD:MATLAB6p5toolboxsymbolicdsolve.matline292x=t-Int

7、(1/(a-a^3+1),a=``..x)+C1=0故只有部分非线性微分方程有解析解。7.2微分方程问题的数值解法7.2.1微分方程问题算法概述微分方程求解的误差与步长问题：7.2.2四阶定步长Runge-Kutta算法及MATLAB实现function[tout,yout]=rk_4(odefile,tspan,y0)％y0初值列向量t0=tspan(1);th=tspan(2);iflength(tspan)<=3,h=tspan(3);%tspan=[t0,th,h]else,h=tspa

8、n(2)-tspan(1);th=tspan(end);end％等间距数组tout=[t0:h:th]';yout=[];fort=tout'k1=h*eval([odefile‘(t,y0)’]);%odefile是一个字符串变量，为表示微分方程的文件名。k2=h*eval([odefile'(t+h/2,y0+0.5*k1)']);k3=h*eval([odefile'(t+h/2,y0+0.5*k2)']);k4=h*eval([odefile'(t+h,y0+k3)']);