第十章常微分方程組求解

常微分方程（組）的符號(hào)求用求常微分方程（組）的通調(diào)用格式一S=dsolve調(diào)用格式二S=dsolve('eqn1','eqn2',用求常微分方程（組）的特調(diào)用格式三S=dsolve線性常微分方程組的解求解齊次線性常微分方程組YAX的主程名為qcxxcwz.m的求解齊次線性常微分方程組YAX的主程序如下function[X,E,V]=qcxxcwz(A)symsx歐拉（Euler）方法的程向前歐 及其誤差估9.2.1

0xyx0h0.0750.0075，并計(jì)算誤差，畫出精確解和數(shù)值解的圖形解編寫并保存名為Eulerli1.m的計(jì)算和畫圖的主程序如functionP=Eulerli1(x0,y0,b,h)n=(b-x0)/h;X=zeros(n,1);Y=zeros(n,1);k=1;X(k)=x0;fork=1:nY(k+1)=Y(k)+h*(X(k)-Y(k));y=X-1+2*exp(-X);plot(X,Y,'mp',X,y,'b-')xlabel('X')ylabel('確解y=x-1+2exp(-x)')解y=x-1+2exp(-x)')>>h1=0.0075;holdholdy=x-1+2exp(-向前歐拉方法的兩種程(一)向前歐 的主程序用向前歐 求解常微分方程初值問題的數(shù)值解及其截?cái)嗾`差的主function[h,k,X,Y,P,REn]=Qeuler1( x=x0;h=(b-x)/n;X=zeros(n,1);y=y0;Y=zeros(n,1);k=1;X(k)=x;Y(k)=y';fork=2:n+1 ifx=x+h;y=y+h*fxy;X(k)=x;Y(k)=y';label('X')ylabel('symsdy2,例9.2.2用向前歐拉求解初值問dy3y8x7,y(0)n10100，并將計(jì)算結(jié)果與精確解作比較，寫出在每個(gè)子區(qū)間[xkxk1解>>subplot(2,1,1) holdonS1=8/3*x1-29/9+38/9*exp(-3*x1),plot(x1,S1,'b-')解holdoffjdwuc1=S1-Y1;jwY1=S1-Y1; holdonS1=8/3*x2-29/9+38/9*exp(-3*x2),plot(x2,S1,'b-')n10100時(shí)，所給的初值問題在[0,1x成的數(shù)組Xi(i=1,2)，利用向前歐拉求出的與Xi(i=1,2)對(duì)應(yīng)的數(shù)值解Yi(i=1,2)，YijwYi(i=1,2)xwYi(i=1,2)（略），每個(gè)子區(qū)間[xkxk1]上的局部截?cái)嗾`差Ren2=1/200*dy2，近似解Y與精確解的圖形.(二)向前歐 的主程序 序2function[k,X,Y,P,REn]=Qeuler2( x=x0;n=fix((b-x)/h);X=zeros(n+1,1);y=y0;Y=zeros(n+1,1);k=1;X(k)=x;Y(k)=y';fork=2:n+1 y=Y(k),k=k+1,plot(X,Y,'rp'),grid,xlabel('X')ylabel('Y')k1=1:n;k=[0,k1];symsdy2,例9.2.3用向前歐 求解初值問題dyy

x,y(0)1n10100將計(jì)算結(jié)果與精確解作比較，寫出在每個(gè)子區(qū)間[xk,xk1]上的局部截?cái)嗾`差，畫出數(shù)值解與精確解在區(qū)間[0,2上的圖形.解>>subplot(2,1,1) holdonS1=1/6*(6+12*X+30*exp(2*X)).^(1/2);plot(X,S1,'b-')值解dy/dx=y-x/(3y)，y(0)=1在[0,2]上的精確解')holdoffjdwucY=S1-Y;jwY=S1-Y; holdondy/dx=y-x/(3y)，y(0)=1在[0,2]上的精確解')運(yùn)行后屏幕顯示取n10,100時(shí)，此初值問題在[0,2上的自變量x的向量X，X1，利用向前歐拉求出的與X，X1對(duì)應(yīng)的數(shù)值解Y，Y1及其絕對(duì)誤差向量jwY，jwY1和相對(duì)誤差向量xwY，xwY1（略），每個(gè)子區(qū)間[xk,xk1]上的局部截?cái)嗾`差Ren2=1/200*dy2，數(shù)Y與精確解的圖形.向后歐拉方法的程用向后歐拉方法求解常微分方程初值問題的數(shù)值解的主程序function[X,Y,n,P]=Heuler1( n=fix((b-x0)/h);X=zeros(n+1,1);k=1;X(k)=x0;Y(k,:)=y0;%繪圖.%產(chǎn)生初值.fori=2:n+1X(i)=x0+h; 主循環(huán)whileWu>tol x0=x0+h;Y(i,:)=y0;plot(X,Y,'ro')gridonxlabel('X')ylabel('值解X=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1;例9.2.4用向后歐拉公式求解區(qū)間[0,2上的初值問題dy3y8xy(0)1h0.05，并與精確解作比較，在同一個(gè)坐標(biāo)系中作出圖形.h0.01h與誤差的關(guān)系. >> holdonS2=8/3*X1-29/9+38/9*exp(-3*X1),plot(X1,S2,'b-holdoff holdonS1=8/3*X-29/9+38/9*exp(-3*X),holdoffS1H,YxiwYjuwY(略).改進(jìn)的歐拉方法的程梯形的程(一)梯 的程 求解常微分方程初值問題的數(shù)值解的主程序1function[X,Y,n,P]=odtixing1( k=1;X(k)=x0;Y(k,:)=y0;%繪圖.%產(chǎn)生初值.fori=2:n+1 主循環(huán)whileWu>tol Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:),Y(i,:)=p1;x0=x0+h;Y(i,:)=y0;plot(X,Y,'ro')gridonxlabel('X')ylabel('值解

X=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1;P=[n',X,Y]例9.3.1用梯 求解區(qū)間[0,2]上的初值問題dy3y8x7,y(0)1，h0.05，精度為101，并與精確解作比較，在同一個(gè)坐標(biāo)系中畫出圖形解>>x0=0;y0=1;b=2;tol=0.1; holdonplot(X,S1,'b-'),holdoffjuwYt=S1-Yt;xiwYt=juwYt./Yt;[0,2上的自變量X處數(shù)值解Yi(i=t,q)和精確解S1H,Yi的相對(duì)誤差xiwYi和絕對(duì)juwYi(略)及其數(shù)值解和精確解的圖形.例10.4.2分別用自適應(yīng)梯形和向前歐拉分別求解區(qū)間[0,2]上的初值問dy3y8x7,y(01，取精度為101作出圖形解>>b=2;tol=1.e-[Ht,X,Yt,k,h,Pt]=odtixing2( holdonS1=8/3*X-29/9+38/9*exp(-3*X),holdoffhold holdoff用自適應(yīng)梯形計(jì)算y，在上的數(shù)值解=，01在上的精確解用向前歐拉計(jì)算=，在上的數(shù)值解值解運(yùn)行后屏幕顯示取精度為101，分別用梯形和向前歐拉求解此初值問題在區(qū)間[0,2XYi(i=t,qS1H等(略).改進(jìn)的歐 的程用自適應(yīng)改進(jìn)的歐 求解常微分方程初值問題的數(shù)值解的主程function ifnargin<5|isempty(tol),tol=1.e-6;ifnargin<6|isempty(trace),trace=0;x=x0;h=0. y=y0(:);p=128;n=fix((b-x0)/h);H=zeros(p,1);X=zeros(p,1);Y=zeros(p,length(y));k=1;X(k)=x;Y(k,:)=y';繪圖主循環(huán)while(x<b)&(x+h>x)ifx+h>b ,x,y);當(dāng)誤差可接受時(shí)重寫解.ifx=x+h;y1=y+h*fxy;fxy1=feval( fxy=fxy(:);y2=(fxy+fxy1)/2;y=y+h*y2;k=k+1;ifk>length(X)X=[X;zeros(p,1)];Y=[Y;zeros(p,length(y))];H(k)=h;X(k)=x;Y(k,:)=y'plot(X,Y,'mh')gridxlabel('自變量X'),ylabel('因變量Y')ifdelta~=0.0ifdisp('Singularitylikely.'),xH=H(1:k);X=X(1:k);Y=Y(1:k,:);n=1:k;例9.4.2分別用梯形公式和改進(jìn)的歐拉公式求解區(qū)間[0,2上的初值問題dy3y8x7,y(01，取精度為101,形解>>x0=0;y0=1;b=2;tol=1.e- holdonS1=8/3*X-29/9+38/9*exp(-3*X),holdoffholdon holdoff','dy/dx=8x-3y-7，y(0)=1在[0,2]上的精確解','用改進(jìn)的歐拉計(jì)算,juwY=S1-Y;xiwY=juwY./Y;運(yùn)行后屏幕顯示取精度為101，分別用梯形和改進(jìn)的歐拉求此解初值問題在區(qū)間[0,2XYtYS1HY和絕對(duì)誤差juwY龍格-庫塔方二階龍格-庫塔方法的程用二階龍格—庫塔方法求解常微分方程初值問題的數(shù)值解的主程 x=x0;y=y0;p=128;k=1;X(k)=x;Y(k,:)=y';%繪圖.fork=2:n+1 ,x,y); grid,xlabel('X'),ylabel('fork=2:n+19.5.1dy1

2xy1

,0x的數(shù)值解，取c114c23

yx0a2b2123h14畫出精確解和數(shù)值解的圖形解在工作窗口輸入下面的程>> ④將運(yùn)行后計(jì)算的結(jié)果列入（略），三階龍格-庫塔方法及其程用三階龍格-庫塔方法求解常微分方程初值問題的數(shù)值解的主程 x=x0;y=y0;p=128;k=1;X(k)=x;Y(k,:)=y';%繪圖.fork=2:n+1 X(k)=x;Y(k,:)=y;grid,xlabel('X')ylabel('fork=2:n+1例9.5.2用常用的三階龍格-庫塔求初值問dy1

2xy1

,0xyx0c116c246c31

并計(jì)算與精確解的誤差，畫出精確解和數(shù)值解的圖形解在工作窗口輸入下面的程>>x0=0;b=2;c1=1/6;C=[c1,c2,c3,a2,a3,b21,b31,b32]; 將運(yùn)行后計(jì)算的結(jié)果（略），四階龍格-庫塔方法及其程用一般四階龍格-庫塔方法式求解常微分方程初值問題的數(shù)值解的主程序 x=x0;y=y0;p=128;n=fix((b-x0)/h);fxy=zeros(p,1);wucha=zeros(p,1);wch=zeros(p,1);X=zeros(p,1);Y=zeros(p,length(y));k=1;X(k)=x;Y(k,:)=y';%繪圖.fork=2:n+1x=x+h;a2=C(5);a4=C(7);b21=C(8);b31=C(9);b32=C(10);b41=C(11);b42=C(12);b43=C(13);c1=C(1);c2=C(2);c3=C(3);x1=x+a2*h;x2=x+a3*h; X(k)=x;Y(k,:)=y;k=k+1plot(X,Y,'rp',X,fxy,'bo')gridxlabel('自變量X'),ylabel('因變量Y')fork=2:n+1X=X(1:k);用常用四階龍格-庫塔方法求解常微分方程初值問題的數(shù)值解的主程序function[x,y]=RKc4(x0,b,y0,h)forx=x0:h:bx=x0+h;[x,y],x0=x;y0=y;例9.5.3用常用的四階龍格-庫塔求初值問dy1

2xy1x2

0xyx0 方法1在工作窗口輸入下面的程>>x0=0;b=2; C=[c1,c2,c3,a2,a3,b21,b31,b32]; c1=1/6;c2=2/6;c3=2/6;c4=1/6;a2=1/2;a3=1/2;C=[c1,c2,c3,c4,a2,a3, 將運(yùn)行后計(jì)算的結(jié)果（略），畫出精確解和數(shù)值解的圖形方法 >>x0=0;b=2;運(yùn)行后輸出結(jié)果（略線性多步法及其程亞當(dāng)斯（Adams）顯式及其程function[k,X,Y,wucha,P]=Adams4x( x=x0;y=y0;p=128;n=fix((b-x0)/h);ifn<5,return,end;Y=zeros(p,length(y));f=zeros(p,1);k=1;X(k)=x;Y(k,:)=y';fork=2:4c4=1/6;a2=1/2;a3=1/2;b32=1/2;b41=0;b42=0;b43=1;x1=x+a2*h;x3=x+a4*h;k1=feval( y1=y+b21*h*k1;x=x+h; y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4);X(k)=x;Y(k,:)=y; for Y(k+1)=Y(k)+(h/24)*((f(k-3:k))'*[-937-59f(k+1)=feval( ,X(k+1),Y(k+1));f(k)=f(k+1);k=k+1;forwucha(k)=norm(Y(k)-Y(k-1));k=k+1;X=X(1:n+1);Y=Y(1:n+1,:);wucha=wucha(1:n,:);例9.6.1先利用常用的四階龍格-庫塔求初值問dy

2xy1x

,0xyx0解在工作窗口輸入下面的程>>y=dsolve('Dy=1->> y=(X+1/3*X.^3)./(1+X.^2);b31=0;b41=0;b42=0;b43=1;a3=1/2;a4=1;b21=1/2;C=[c1,c2,c3,c4,a2,a3, grid,xlabel('X')ylabel('Y')wchY=abs(y-Y),wchY1=abs(y-Y1),forn=1:k,m(1,n)=n-1,end,運(yùn)行后屏幕顯示圖和計(jì)算結(jié)果（略亞當(dāng)斯隱式及其程例9.6.2先利用常用的四階龍格-庫塔求初值問dyxy,0xyx1ex的誤差，在同一圖形窗口畫出精確解和數(shù)值解的圖形.解編寫并保存名為Adams4y.m的程function[k,X,Y,wucha,P]=Adams4y(x0,b,y0,h)x=x0;y=y0;p=128;ifn<5,Y=zeros(p,length(y));f=zeros(p,1);k=1;X(k)=x;Y(k,:)=y';forx1=x+h/2;x2=x+h/2;X(k)=x;Y(k,:)=y;k=k+1;fork=3:nfor>>[k,X,Y,wucha,P]=Adams4y(x0,b,y0,h)b31=0;b41=0;a2=1/2;a3=1/2;C=[c1,c2,c3,c4,a2,a3, xlabel('X')ylabel('','dy/dx=x-y，y(0)=0的精確解y=x-1+exp(-x)')wchY1=abs(y-Y1),m=zeros(1,k),forn=1:k運(yùn)行后屏幕顯示圖和計(jì)算結(jié)果（略（二）改進(jìn)的亞當(dāng)斯方法及其程用改進(jìn)的亞當(dāng)斯方法求解常微分方程初值問題的數(shù)值解的主程序function[k,X,Y,wucha,P]=Adams4y1( x=x0;y=y0;p=128;n=fix((b-x0)/h);ifn<5,X=zeros(p,1);Y=zeros(p,length(y));f=zeros(p,1);k=1;X(k)=x;Y(k,:)=y';fork=2:4x2=x+a3*h;x3=x+a4*h; fork=4:n Y(k+1)=Y(k)+(h/24)*((f(k-2:k+1))'*[1-5199]');f(k+1)=feval( f(k)=f(k+1);k=k+1;for例 用改進(jìn)的亞當(dāng)斯方法和常用的四階龍格-庫塔求例9.6.2中初值問解