用MATLAB求解微分方程及微分方程組.ppt

1.微分方程的解析解,求微分方程（組）的解析解命令:,dsolve(方程1,方程2,方程n,初始條件,自變量),運(yùn)行結(jié)果：u=tan(t-c),用MATLAB求解微分方程,解輸入命令：dsolve(Du=1+u2,t),解輸入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),運(yùn)行結(jié)果為:y=3e-2xsin（5x）,解輸入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%將x化簡(jiǎn)y=simple(y)z=simple(z),運(yùn)行結(jié)果為：x=(c1-c2+c3+c2e-3t-c3e-3t)e2ty=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t,2.用Matlab求常微分方程的數(shù)值解,t，x=solver（f,ts,x0,options）,1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫成.,2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.,注意:,解:令y1=x，y2=y1,1、建立m-文件vdp1000.m如下：functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);,2、取t0=0，tf=3000，輸入命令：T,Y=ode15s(vdp1000,03000,20);plot(T,Y(:,1),-),3、結(jié)果如圖,解1、建立m-文件rigid.m如下：functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,2、取t0=0，tf=12，輸入命令：T,Y=ode45(rigid,012,011);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),3、結(jié)果如圖,圖中，y1的圖形為實(shí)線，y2的圖形為“*”線，y3的圖形為“+”線.,導(dǎo)彈追蹤問(wèn)題,設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于x軸上點(diǎn)A(1,0)處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦.如果乙艦以最大的速度v0(是常數(shù))沿平行于y軸的直線行駛，導(dǎo)彈的速度是5v0，求導(dǎo)彈運(yùn)行的曲線方程.又乙艦行駛多遠(yuǎn)時(shí)，導(dǎo)彈將它擊中？,解法一（解析法）,由(1),(2)消去t整理得模型:,解法二(數(shù)值解),1.建立m-文件eq1.mfunctiondy=eq1(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)2)/(1-x);,2.取x0=0，xf=0.9999，建立主程序ff6.m如下：x0=0，xf=0.9999x,y=ode15s(eq1,x0 xf,00);plot(x,y(:,1),b.)holdony=0:0.01:2;plot(1,y,b*),結(jié)論:導(dǎo)彈大致在（1，0.2）處擊中乙艦,令y1=y,y2=y1，將方程（3）化為一階微分方程組。,解法三(建立參數(shù)方程求數(shù)值解),設(shè)時(shí)刻t乙艦的坐標(biāo)為(X(t)，Y(t)，導(dǎo)彈的坐標(biāo)為(x(t)，y(t).,3因乙艦以速度v0沿直線x=1運(yùn)動(dòng)，設(shè)v0=1，則w=5，X=1，Y=t,4.解導(dǎo)彈運(yùn)動(dòng)軌跡的參數(shù)方程,建立m-文件eq2.m如下：functiondy=eq2(t,y)dy=zeros(2,1);dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2);dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2);,取t0=0，tf=2，建立主程序chase2.m如下：t,y=ode45(eq2,02,00);Y=0:0.01:2;plot(1,Y,-)，holdonplot(y(:,1),y(:,2),*),軌跡圖如下,例：飲酒模型,模型1：快速飲酒后，胃中酒精含量的變化率,模型2：快速飲酒后，體液中酒精含量的變化率,即,用Matlab求解模型2：,symsxyk1k2Mtx=dsolve(Dx+k2*x=k1*M*exp(-k1*t),x(0)=0,t),運(yùn)行結(jié)果：,M*k1/(-k1+k2)*exp(-k2*t+t*(-k1+k2)-exp(-k2*t)*M*k1/(-k1+k2),即：,用以下一組數(shù)據(jù)擬合上述模型中的參數(shù)k1、k2：,M=64000/490=130.6122（毫克百毫升）,建立函數(shù)文件：functionf=curvefun1(k,t)f=k(1)*64000/490*(exp(-k(2)*t)-exp(-k(1)*t)/(k(1)-k(2),輸入擬合數(shù)據(jù)：t=00.250.50.7511.522.533.544.55678910111213141516;c=03068758282776868585150413835282518151210774;,任取k1、k2的一組初始值：k0=2,1;,輸入命令：k=lsqcurvefit(curvefun1,k0,t,c),運(yùn)行結(jié)果為：,k=1.32400.2573,作圖表示求解結(jié)果：,t1=0:0.1:18;f=curvefun1(k,t1);plot(t,c,ko,t1,f,r-),模型2：慢速飲酒時(shí)，體液中酒精含量的變化率,則有；,其中,M為飲酒的總量，T為飲酒的時(shí)間,