微分方程組Matlab解析解數(shù)值解

1、求微分方程的數(shù)值解一常微分方程數(shù)值解的定義二建立數(shù)值解法的一些途徑三用Matlab軟件求常微分方程的數(shù)值解返 回微分方程的解析解 求微分方程（組）的解析解命令:dsolve(方程1, 方程2,方程n, 初始條件, 自變量)To Matlabff1 結(jié) 果：u = tg(t-c)解 輸入命令 ： 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化簡 y=simple(y) z=simple(z)結(jié) 果 為：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e

2、-4t+c2e-4t+c2e-3t-c3e-3t+(c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t To Matlabff3返 回微分方程的數(shù)值解一常微分方程數(shù)值解的定義 在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實際上對初值問題，一般是要求得到解在假設(shè)干個點上滿足規(guī)定精確度的近似值，或者得到一個滿足精確度要求的便于計算的表達(dá)式。因此，研究常微分方程的數(shù)值解法是十分必要的。返 回2、使用數(shù)值積分對方程y=f(x,y), 兩邊由xi到xi+1積分，并利用梯形公式，有：實際應(yīng)用時，與歐拉公式結(jié)合使用：此即改進(jìn)的歐拉法。故有公式：3、

3、使用泰勒公式 以此方法為根底，有龍格-庫塔法、線性多步法等方法。4、數(shù)值公式的精度 當(dāng)一個數(shù)值公式的截斷誤差可表示為Ohk+1時k為正整數(shù)，h為步長，稱它是一個k階公式。k越大，那么數(shù)值公式的精度越高。歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。返 回 1、在解n個未知函數(shù)的方程組時，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫成. 2、使用Matlab軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組.注意:解: 令 y1=x，y2=y11、建立m-文件vdp1000.m如下： functi

4、on dy=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,0 3000,2 0); plot(T,Y(:,1),-)3、結(jié)果如圖To Matlabff4導(dǎo)彈追蹤問題 設(shè)位于坐標(biāo)原點的甲艦向位于x軸上點A(1, 0)處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對準(zhǔn)乙艦.如果乙艦以最大的速度v0(是常數(shù))沿平行于y軸的直線行駛，導(dǎo)彈的速度是5v0，求導(dǎo)彈運行的曲線方程.又乙艦行駛多遠(yuǎn)時，導(dǎo)彈將它擊中？解法一解析法解法二(數(shù)值解)1.建立

5、m-文件eq1.m function dy=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.9999 x,y=ode15s(eq1,x0 xf,0 0); plot(x,y(:,1),b.) hold on y=0:0.01:2; plot(1,y,b*) 結(jié)論: 導(dǎo)彈大致在1，0.2處擊中乙艦To Matlab(ff6)令y1=y,y2=y1，將方程3化為一階微分方程組。解法三(建立參數(shù)方程求數(shù)值解) 設(shè)時刻t乙艦的坐標(biāo)為

6、(X(t)，Y(t)，導(dǎo)彈的坐標(biāo)為(x(t)，y(t).3因乙艦以速度v0沿直線x=1運動，設(shè)v0=1，那么w=5，X=1，Y=t5. 結(jié)果見圖1導(dǎo)彈大致在1，0.2處擊中乙艦，與前面的結(jié)論一致.圖1圖2返 回 在chase2.m中，按二分法逐步修改tf，即分別取tf=1，0.5,0.25,直到tf=0.21時，得圖2.結(jié)論：時刻t=0.21時，導(dǎo)彈在1，0.21處擊中乙艦。To Matlab(chase2)2. 模型求解(1) w=20時,建立m-文件eq3.m如下: function dy=eq3(t,y) dy=zeros(2,1); dy(1)=20*(10+20*cos(t)-y(1

7、)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); dy(2)=20*(20+15*sin(t)-y(2)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);取t0=0，tf=10，建立主程序chase3.m如下： t0=0;tf=10; t,y=ode45(eq3,t0 tf,0 0); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,-) hold on plot(y(:,1),y(:,2),*) 在chase3.m，不斷修改tf的值,分別取

8、tf=5, 2.5, 3.5,至3.15時,狗剛好追上慢跑者.To Matlab(chase3)地中海鯊魚問題 意大利生物學(xué)家Ancona曾致力于魚類種群相互制約關(guān)系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚等的比例有明顯增加見下表，而供其捕食的食用魚的百分比卻明顯下降.顯然戰(zhàn)爭使捕魚量下降，食用魚增加，鯊魚等也隨之增加，但為何鯊魚的比例大幅增加呢？ 他無法解釋這個現(xiàn)象，于是求助于著名的意大利數(shù)學(xué)家V.Volterra，希望建立一個食餌捕食系統(tǒng)的數(shù)學(xué)模型，定量地答復(fù)這個問題.首先，建立m-文件shier.m如下： function dx=shier(t,x) dx=zeros(2,1); dx(1)=x(1)*(1-0.1*x(2); dx(2)=x(2)*(-0.5+0.02*x(1);其次，建立主程序shark.m如下： t,x=ode45(shier,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2)To Matlab(shark)模型求解:1、分別用m-文件shier1.m和sh