第七章微分方程问题的解法 微分方程的解析解方法 ·微分方程问题的数值解法 微分方程问题算法概述 四阶定步长 Runge- Kutta算法及 MATLAB实现 一阶微分方程组的数值解 微分方程转换 特殊徼分方程的数值解
第七章 微分方程问题的解法 • 微分方程的解析解方法 • 微分方程问题的数值解法 –微分方程问题算法概述 –四阶定步长 Runge-Kutta算法及 MATLAB 实现 –一阶微分方程组的数值解 –微分方程转换 • 特殊微分方程的数值解
7.1微分方程的解析解方法 格式: y-dsolve(fi, f 格式:指明自变量 y=dsolve(f, f2 5…··m f即可以描述微分方程,又可描述初始条件 或边界条件。如: 描述微分方程时y4)()=7→D4y=7 描述条件时 j(2)=3→D2y(2)=3
7.1 微分方程的解析解方法 • 格式： y=dsolve(f1 , f2 , …, fm) • 格式：指明自变量 y=dsolve(f1 , f2 , …, fm ,’x’) fi即可以描述微分方程，又可描述初始条件 或边界条件。如： 描述微分方程时 描述条件时 (4) y t D y ( ) 7 4 7 = = y D y (2) 3 2 (2) 3 = =
例:输入信号为u(1)=e-51cos(2+1)+5, 求该微分方程的通解 4)(1)+10y()+35y"(+50y′()+24y()=5(+4l(0)+2( step > syms t; u=exp(-5*t)*coS(2*t+1)+5 >>uu=5*diff(u, t, 2)+4 diff(u, t)+2*u uu 87*exp(-5*1)*cO(2*计1)+92*exp(-5*t)*sin(2*t廿+1)+10 step syms t y; >>y= dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y-12 87*exp(-5*1)*cOS(2*t+1)+92*Xp(-5*1)*sin(2*t+1)+10])
例： >> syms t; u=exp(-5*t)*cos(2*t+1)+5; >> uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu = 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10 >> syms t y; >> 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().5343-5c0s(2t+1)520s1n(2t+1)+ 5475r 12520 CIe+C 2e-3t +c 3e-2 + cae 假设已知y(0)=3,y(0)=2,y(O)=y"(O)=0 >y=dsolve(lD4y+10 D3y+35*D2y+50*Dy+24 y= 87*eXp(-5*1)*cOS(2*+1)+92*Xp(-5*)*n(2*t+1) +10,y(0)=3,"Dy(0)=2,"D2y(0=0,D3y(0=0 5343-5tcos(21+1 547 e )-e-5sin(2t+1)+ 2520 520 445 51 27 41 cos 1 SIn sin 1 26 13 30 COS 1+ 179 cos 1+-sin 1+ 133 8 8 3 30cos1≠ sin 1+19e 60
>> 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)=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%ratO最接近有理数的分数 方程的解: 5343-5 coS(2t +1) 5475t e sin(2t+1) 12520 520 8704-2243-4tx136e 5025 2981 185 254 131 判断误差: vpa(-445/26*cos(sym(1)-51/13*sin(1)-692+8704/185) ans 114731975864790922564144636e-4
分别处理系数，如： >> [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))] ans = -8704 185 % rat()最接近有理数的分数 判断误差： >> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185) ans = .114731975864790922564144636e-4
设置y(0)=1/2,y(x)=1,y"(2x)=0,y'(2x)=1/5: y-dsolve(ld4y10*D3y+35 D2y+50*Dy+24*y= 87*exp(-5*1)*COS(2*t+1)+92*eXp(-5*1)*sin(2*t+1)+ 10]y(0)=12,Dy(pi)=1,D2y(2*pi)=0,Dy(2*pi)=15) 如果用推导的方法求Ci的值,每个系数的解析解至少要写出10 数行,故可采用有理式近似的方式表示 >>vpa(y10)%有理近似值 ans 196361839*eXp(-5*t)+.416666667 4785447354*sin(t)*cos(t*exp(-5*t)-.4519262218e 1*cos(2.*t)*exp(-5.*t2.392723677*cos(t)^2*exp( 5*t)+2259631109sin(2.*t)*exp(-5.*t473690.0893*eXp( 3.*t)+3131963786*eXp(-2.*t)-219.1293619eXp 1.*t)+442590.9059eXp(-4.*t)
>> 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.196361839*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.*t)
例:求x()=x()(1-x2(m)的解析解 s> syms tX X-dsolve(DX- [1/(1+exp(-2*t)*C1)^(1/2) [-1/(1+exp(-2*t)*C1)(12) 改变原微分方程:x()=x((1-x2(m)+1 syms tx; x-dsolve(Dx=x*(1-x2)+1) Warning: Explicit solution could not be found; implicit solution returned > In D: MATLAB6p5 toolbox symbolicldsolve m at line 292 Int(1a^3+1),a=3.xy)+Cl=0 故只有部分非线性微分方程有解析解
• 例： >> syms t x >> x=dsolve('Dx=x*(1-x^2)') x = [ 1/(1+exp(-2*t)*C1)^(1/2)] [ -1/(1+exp(-2*t)*C1)^(1/2)] >> syms t x; x=dsolve('Dx=x*(1-x^2)+1') Warning: Explicit solution could not be found; implicit solution returned. > In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292 x = t-Int(1/(a-a^3+1),a=``..x)+C1=0 故只有部分非线性微分方程有解析解
7.2微分方程问题的数值解法 7.2.1微分方程问题算法概述 (t)=f(t,m() 其中x(1)=[x1(,x2(),…,xn()为状态向量 f()=[f1(,1(),…,()可以是任意非线性函数
7.2 微分方程问题的数值解法 7.2.1 微分方程问题算法概述
微分方程求解的误差与步长问题: (t)=f(t,(t) Euer算法 0时刻系统状态向量的值为(t0) 微分方程左侧的导数近似为: (ae(to +h)-a(to))/(to +h-to (0+h时刻微分方程的近似解为 e(to +h)= c(to)+hf(to, a(to)
微分方程求解的误差与步长问题：
在t+h时刻系统状态向量的值为 a(to +h)=a(to +h)+ ro= o +hf(t, co)+ Ro 简记x1=(0+h 在tk时刻系统的状态向量为k 在tk+h时刻 Euler算法的数值解为: Ck+1=k +hf(t, k