常微分方程組求解

1、第十章常微分方程(組)求解J10.1 常微分方程(組)的MATLA箭號求解10.1.1 用MATLA歌常微分方程(組)的通解調(diào)用格式一：S=dsolve('eqn','var')調(diào)用格式二：S=dsolve('eqn1','eqn2',.,'eqnm','var')10.1.2 用MATLA求常微分方程(組)的特解調(diào)用格式三：S=dsolve('eqn'，匕ondition1',，'conditionn','var')調(diào)用格式四：S=dsolv

2、e('eqn1','eqn2',.,'eqnm','condition1','condition2''var')10.1.3 線性常微分方程組的解法求解齊次線性常微分方程組Y=AX的MATLA在程序名為qcxxcwz.m的求解齊次線性常微分方程組Y=AXMATLAB主程序如下:functionX,E,V=qcxxcwz(A)symsxE=eig(A);V,n=eig(A);X=exp(E*x)'*V;10.3歐拉(Euler)方法的MATLAB?序向前歐拉公式及其誤差估計例10.2.1用歐拉

3、方法求初值問題idy二x-y,0-x-1,dxiyXz0=1的數(shù)值解，分別取h=0.0750,0.0075,并計算誤差，畫出精確解和數(shù)值解的圖形解編寫并保存名為Eulerli1.m的MATLAB十算和畫圖的主程序如下functionP=Eulerli1(x0,y0,b,h)n=(b-x0)/h;X=zeros(n,1);Y=zeros(n,1);k=1;X(k)=x0;Y(k)=y0;fork=1:nX(k+1)=X(k)+h;丫(k+1)=丫(k)+h*(X(k)-Y(k);k=k+1;endy=X-1+2*exp(-X);plot(X,Y,'mp',X,y,'b-&

4、#39;)gridxlabel('自變量X'),ylabel('因變量Y')title(解y=x-1+2exp(-x)'legend(解y=x-1+2exp(-x)''用向前歐拉公式求dy/dx=x-y,y(0)=1在0,1上的數(shù)值解和精確)'h=0.075時，dy/dx=x-y,y(0)=1在0,1上的數(shù)值解','精確)jwY=y-Y;xwY=jwY./y;k1=1:n;k=0,k1;P=k',X,Y,y,jwY,xwY;在MATLAB：作窗口輸入下面的程序>>x0=0;y0=1;b=1;h=

5、0.0750;P=Eulerli1(x0,y0,b,h)在MATLAB：作窗口輸入下面的程序>>h1=0.0075;P1=Eulerli1(x0,y0,b,h1)legend('h1=0.0075時，dy/dx=x-y,y(0)=1在0,1上的數(shù)值解','精確解y=x-1+2exp(-x)')向前歐拉方法的三種MATLA毓序(一)向前歐拉公式的MATLABE程序1用向前歐拉公式(10.8)求解常微分方程初值問題(10.5)的數(shù)值解及其截斷誤差公式的MATLA在程序1functionh,k,X,Y,P,REn=Qeuler1(funfcn,x0,y0,

6、b,n,tol)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+1fxy=feval(funfcn,x,y);delta=norm(h*fxy,'inf);wucha=tol*max(norm(y,'inf),1.0);ifdelta>=wuchax=x+h;y=y+h*fxy;X(k)=x;Y(k)=y'endplot(X,Y,'rp')gridlabel('自變量X'),ylabel('因變量Y')title

7、('用向前歐拉(Euler)公式計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解')endPX,Y;symsdy2,REn=0.5*dy2*hA2;例10.3.2用向前歐拉公式(10.8)求解初值問題dy=-3y+8x-7,y(0)=1，dx分別取n=10,100,并將計算結(jié)果與精確解作比較，寫出在每個子區(qū)間xk,xk+上的局部截斷誤差公式，畫出數(shù)值解與精確解在區(qū)間0,1上的圖形.解輸入程序>>subplot(2,1,1)x0=0;y0=1;b=1-1.e-4;n=100;tol=1.e-4;h1,k1,x1,Y1,P1,Ren1=QEuler1(

8、unfcn,x0,y0,b,n,tol)holdonS1=8/3*x1-29/9+38/9*exp(-3*x1),plot(x1,S1,'b-')title('用向前歐拉公式計算dy/dx=8x-3y-7,y(0)=1在0,1上的數(shù)值解')legend('n=100時，dy/dx=8x-3y-7,y(0)=1在0,1上的數(shù)值解','dy/dx=8x-3y-7,y(0)=1在0,1上的精確解')holdoffjdwuc1=S1-Y1;jwY1=S1-Y1;xwY1=jwY1./S1;k1=1:n;k=0,k1;P1=k',x

9、1,Y1,S1,jwY1,xwY1subplot(2,1,2)n1=10;h2,k2,x2,Y2,P2,Ren2=QEuler1(unfcn,x0,y0,b,n1,tol)holdonS1=8/3*x2-29/9+38/9*exp(-3*x2),plot(x2,S1,'b-')legend('n=10時，dy/dx=8x-3y-7,y(0)=1在0,1上的數(shù)值解','dy/dx=8x-3y-7,y(0)=1在0,1上的精確解')holdoffjwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=0,k1;P2=k',x2,

10、Y2,S1,jwY2,xwY2運行后屏幕顯示分別取n=10,100時，所給的初值問題在0,1上的自變量x的值構(gòu)成的數(shù)組Xi(i=1,2),利用向前歐拉公式(10.8)求出的與Xi(i=1,2)對應(yīng)的數(shù)值解Yi(i=1,2),Yi的絕對誤差jwYi(i=1,2)和相對誤差xwYi(i=1,2)(略)，每個子區(qū)間xxk書上的局部截斷誤差公式Ren2=1/200*dy2,近似解丫與精確解的圖形.(二)向前歐拉公式的MATLA更程序2用向前歐拉公式(10.8)求解常微分方程初值問題(10.5)的數(shù)值解及其截斷誤差公式的MATLA在程序2functionk,X,Y,P,REn=Qeuler2(funfc

11、n,x0,y0,b,h)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+1X(k)=x+(k-1)*h,fxy=feval(funfcn,x,y),Y(k)=y+h*fxyy=Y(k),k=k+1,plot(X,Y,'rp'),grid,xlabel('自變量X'),ylabel('因變量Y')title('用向前歐拉(Euler)公式計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解')end

12、k1=1:n;k=0,k1;P=k',X,Y;symsdy2,REn=0.5*dy2*hA2;例10.3.4用向前歐拉公式(10.8)求解初值問題一=y,y(0)=1,取dx3yn=10,100,并將計算結(jié)果與精確解作比較，寫出在每個子區(qū)間xk,xk由上的局部截斷誤差公式，畫出數(shù)值解與精確解在區(qū)間0,2上的圖形.解輸入程序>>subplot(2,1,1)x0=0;y0=1;b=2;n=10;h=2/10;k,X,Y,P,REn=Qeuler2(unfcn,x0,y0,b,h)holdonS1=1/6*(6+12*X+30*exp(2*X)A(1/2);plot(X,S1,&

13、#39;b-')title('用向前歐拉公式計算dy/dx=y-x/(3y),y(0)=1在0,2上的數(shù)值解)legend('n=10時，dy/dx=y-x/(3y),y(0)=1在0,2上的數(shù)值解','dy/dx=y-x/(3y),y(0)=1在0,2上的精確解')holdoffjdwucY=S1-Y;jwY=S1-Y;xwY=jwY./Y;k1=1:n;k=0,k1;P1=k',X,Y,S1,jwY,xwYsubplot(2,1,2)n1=100;h1=2/100;k,X1,Y1,P1,Ren1=Qeuler2(unfcn,x0,y0

14、,b,h1)holdonS2=1/6*(6+12*X1+30*exp(2*X1)(1/2);plot(X1,S2,'b-')legend('n=100時，dy/dx=y-x/(3y),y(0)=1在0,2上的數(shù)值解','dy/dx=y-x/(3y),y(0)=1在0,2上的精確解')holdoffjwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:n1;k=0,k1;P2=k',X1,Y1,S2,jwY1,xwY1運行后屏幕顯示取n=10,100時，此初值問題在0,2上的自變量x的向量X,X,利用向前歐拉公式(10.8)求出的與X

15、,X對應(yīng)的數(shù)值解Y,Y及其絕對誤差向量jwY,jwY1和相對誤差向量xwYxwY(略)，每個子區(qū)間xk,xk書上的局部截斷誤差公式Ren2=1/200*dy2,數(shù)值解丫與精確解的圖形.(三)自適應(yīng)向前歐拉公式的MATLA更程序用自適應(yīng)向前歐拉公式求解常微分方程初值問題(10.5)的數(shù)值解的MATLA駐程序2functionH,X,Y,k,h,P=QEuler(funfcn,x0,b,y0,tol)%始化.pow=1/3;ifnargin<5|isempty(tol),tol=1.e-6;end;ifnargin<6|isempty(trace),trace=0;end;x=x0;h

16、=0.0078125*(b-x);y=y0(:);p=128;H=zeros(p,1);X=zeros(p,1);Y=zeros(p,length(y);k=1;X(k)=x;Y(k,:)=y'%繪圖.clc,x,h,yend%主循環(huán).while(x<b)&(x+h>x)ifx+h>bh=b-x;end%計算斜率.fxy=feval(funfcn,x,y);fxy=fxy(:);%計算誤差，設(shè)定可接受誤差.delta=norm(h*fxy,'inf);wucha=tol*max(norm(y,'inf),1.0);%當(dāng)誤差可接受時重寫解.ifd

17、elta<=wuchax=x+h;y=y+h*fxy;k=k+1;ifk>length(X)X=X;zeros(p,1);Y=Y;zeros(p,length(y);H=H;zeros(p,1);endH(k)=h;X(k)=x;Y(k,:)=y'plot(X,Y,'rp'),gridxlabel('自變量X'),ylabel('因變量Y')title('用向前歐拉(Euler)公式計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解')end%更新步長.ifdelta=0.0h=min(h*8,

18、0.9*h*(wucha/deltaFpow);endendif(x<b)disp('Singularitylikely.'),xendH=H(1:k);X=X(1:k);Y=Y(1:k,:);n=1:k;P=n',H,X,Y10.3.4向后歐拉方法的MATLAB?序用向后歐拉方法求解常微分方程初值問題(10.5)的數(shù)值解的MATLA在程序functionX,Y,n,P=Heuler1(funfcn,x0,b,y0,h,tol)n=fix(b-x0)/h);X=zeros(n+1,1);Y=zeros(n+1,1);k=1;X(k)=x0;Y(k,:)=y0;Y1

19、(k,:)=y0;%繪圖.clc,x0,h,y0%產(chǎn)生初值.fori=2:n+1X(i)=x0+h;Y(i,:)=y0+h*feval(funfcn,x0,y0);Y1(i,:)=y0+h*feval(funfcn,X(i),Y(i,:);%主循環(huán).Wu=abs(Y1(i,:)-Y(i,:);whileWu>tolp=Y1(i,:);Y1(i,:)=y0+h*feval(funfcn,X(i),p);Y(i,:)=p;endx0=x0+h;y0=Y1(i,:);Y(i,:)=y0;plot(X,Y,'ro')gridonxlabel('自變量X'),yla

20、bel('因變量Y')title('用向后歐拉公式計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解)endX=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1;P=n',X,Y例10.3.6用向后歐拉公式求解區(qū)間0,2上的初值問題y=3y+8x7,dxy(0)=1的數(shù)值解，取步長h=0.05,并與精確解作比較，在同一個坐標(biāo)系中作出圖形然后再取h=0.01,觀察數(shù)值解與精確解誤差的變化，說明h與誤差的關(guān)系.解輸入程序>>S1=dsolve('Dy=8*x-3*y-7','y(0)=1',

21、9;x')>>x0=0;y0=1;b=2;tol=1.e-1;subplot(2,1,1)h1=0.01;X1,Y1,n,P1=Heuler1(unfcn,x0,b,y0,h1,tol)holdonS2=8/3*X1-29/9+38/9*exp(-3*X1),plot(X1,S2,'b-')legend('h=0.01用向后歐拉公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解','dy/dx=8x-3y-7,y(0)=1在0,2上的精確解')holdoffjuwY1=S2-Y1;xiwY1=juwY1./Y1

22、;L=P1,S2,juwY1,xiwY1subplot(2,1,2)h=0.05;X,Y,n,P=Heuler1(unfcn,x0,b,y0,h,tol)holdonS1=8/3*X-29/9+38/9*exp(-3*X),plot(X,S1,'b-')legend('h=0.05用向后歐拉公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解','dy/dx=8x-3y-7,y(0)=1在0,2上的精確解')holdoffjuwY=S1-Y;xiwY=juwY./Y;L=P,S1,juwY,xiwY運行后屏幕顯示用向后歐拉公式計算此

23、初值問題在0,2上的自變量X處數(shù)值解Y和精確解Si及其圖形，步長HY的相對誤差xiwY和絕對誤差juwY(略).10.4改進的歐拉方法的MATLABS序10.4.2梯形公式的兩種MATLA醒序(一)梯形公式的MATLA典序用梯形公式求解常微分方程初值問題(10.5)的數(shù)值解的MATLA駐程序1functionX,Y,n,P=odtixing1(funfcn,x0,b,y0,h,tol)n=fix(b-x0)/h);X=zeros(n+1,1);Y=zeros(n+1,1);k=1;X(k)=x0;Y(k,:)=y0;Y1(k,:)=y0;%繪圖.clc,x0,h,y0%產(chǎn)生初值.fori=2:

24、n+1X(i)=x0+h;fx0y0=feval(funfcn,x0,y0);Y(i,:)=y0+h*fx0y0;fxiyi=feval(funfcn,X(i),Y(i,:);丫1(i,:)=y0+h*(fxiyi+fx0y0)/2;%主循環(huán).Wu=abs(Y1(i,:)-Y(i,:);whileWu>tolP=Y1(i,：),fxip=feval(funfcn,X(i),p);Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:),Y(i,:)=p1;endx0=x0+h;y0=Y1(i,:);Y(i,:)=y0;plot(X,Y,'ro')gri

25、donxlabel('自變量X'),ylabel('因變量Y')title('用梯形公式計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解)endX=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1;P=n',X,Y例10.4.1用梯形公式求解區(qū)間0,2上的初值問題電=_3y+8x_7,y(0)=1,dx取步長h=0.05,精度為10,，并與精確解作比較，在同一個坐標(biāo)系中畫出圖形解輸入程序>>x0=0;y0=1;b=2;tol=0.1;h=0.05;X,Yt,n,Pt=odtixing1(unfcn,x0,b

26、,y0,h,tol)holdonS1=8/3*X-29/9+38/9*exp(-3*X);plot(X,S1,'b-'),holdofflegend('h=0.05,用梯形公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解，dy/dx=8x-3y-7,y(0)=1在0,2上的精確解')juwYt=S1-Yt;xiwYt=juwYt./Yt;Lt=Pt,S1,juwYt,xiwYt運行后屏幕顯示取精度為10，分別用梯形公式和向前歐拉公式求解此初值問題在區(qū)間0,2上的自變量X處數(shù)值解Y(i=t,q)和精確解S,步長H,Y的相對誤差xiwYi和絕對誤差

27、juwY(略)及其數(shù)值解和精確解的圖形.(二)自適應(yīng)梯形公式的MATLAB序用自適應(yīng)梯形公式求解常微分方程初值問題(10.5)的數(shù)值解的MATLA在程序functionH,X,Y,k,h,P=odtixing2(funfcn,x0,b,y0,tol)%始化.pow=1/3;ifnargin<5|isempty(tol),tol=1.e-6;end;ifnargin<6|isempty(trace),trace=0;end;x=x0;h=0.0078125*(b-x);y=y0(:);p=128;n=fix(b-x0)/h);H=zeros(p,1);X=zeros(p,1);Y=z

28、eros(p,length(y);k=1;X(k)=x;Y(k,:)=y'%繪圖.clc,x,h,y%主循環(huán).while(x<b)&(x+h>x)ifx+h>bh=b-x;end%十算余率.fxy=feval(funfcn,x,y);fxy=fxy(:);%計算誤差，設(shè)定可接受誤差.delta=norm(h*fxy,'inf);wucha=tol*max(norm(y,'inf),1.0);%當(dāng)誤差可接受時重寫解.ifdelta<=wuchax=x+h;y1=y+h*fxy;fxy1=feval(funfcn,x,y1);fxy=fxy(

29、:);y2=y+h*fxy1;y=(y1+y2)k=k+1;ifk>length(X)X=X;zeros(p,1);Y=Y;zeros(p,length(y);H=H;zeros(p,1);endH(k)=h;X(k)=x;Y(k,:)=y'plot(X,Y,'go'),gridxlabel('自變量X'),ylabel('因變量Y')title('用自適應(yīng)梯形公式計算dy/dx=f(x,y),y(x0)=y0在x0,b上的end%庾新步長.ifdelta=0.0h=min(h*8,0.9*h*(wucha/delta)Ap

30、ow);endendif(x<b)disp('Singularitylikely.'),xendH=H(1:k);X=X(1:k);Y=Y(1:k,:);n=1:k;P=n',H,X,Y例10.4.2分別用自適應(yīng)梯形公式和向前歐拉公式分別求解區(qū)間0,2上的初值問題曳=-3y+8x-7,y(0)=1,取精度為10,，并與精確解作比較，在同一個坐標(biāo)系中dx作出圖形.解輸入程序>>x0=0;y0=1;b=2;tol=1.e-1;Ht,X,Yt,k,h,Pt=odtixing2(funfcn,x0,b,y0,tol),holdonS1=8/3*X-29/9+3

31、8/9*exp(-3*X),plot(X,S1,'b-'),holdoffholdon,Hq,X,Yq,k,h,Pq=QEuler(unfcn,x0,b,y0,tol),holdoffgridlegend('用自適應(yīng)梯形公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解，dy/dx=8x-3y-7,y(0)=1在0,2上的精確解','用向前歐拉公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解')title('用自適應(yīng)梯形公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解)運行后屏幕顯示取精

32、度為10,，分別用梯形公式和向前歐拉公式求解此初值問題在區(qū)間0,2上的自變量X處數(shù)值解Y(i=t,q)和精確解S及其圖形，步長H等(略).10.4.4改進的歐拉公式的MATLAB?序用自適應(yīng)改進的歐拉公式求解常微分方程初值問題(10.5)的數(shù)值解的MATLA在程functionH,X,Y,k,h,P=gaiEuler(funfcn,x0,b,y0,tol)%始化.pow=1/3;ifnargin<5|isempty(tol),tol=1.e-6;end;ifnargin<6|isempty(trace),trace=0;end;x=x0;h=0.0078125*(b-x);y=y0

33、(:);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'%繪圖.clc,x,h,yend%主循環(huán).while(x<b)&(x+h>x)ifx+h>bh=b-x;end%十算余率.fxy=feval(funfcn,x,y);fxy=fxy(:);%計算誤差，設(shè)定可接受誤差.delta=norm(h*fxy,'inf);wucha=tol*max(norm(y,'inf),1.0);%當(dāng)誤差可接受時重寫解.ifdelta&l

34、t;=wuchax=x+h;y1=y+h*fxy;fxy1=feval(funfcn,x,y1);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=H;zeros(p,1);endH(k)=h;X(k)=x;Y(k,:)=y'plot(X,Y,'mh'),gridxlabel('自變量X'),ylabel('因變量Y')title('用改進的歐拉公式計算dy/dx=f(x,y),y(x0)

35、=y0在x0,b上的數(shù)值解')end%更新步長.ifdelta=0.0h=min(h*8,0.9*h*(wucha/deltaFpow);endendif(x<b)disp('Singularitylikely.'),xendH=H(1:k);X=X(1:k);Y=Y(1:k,:);n=1:k;P=n',H,X,Y例10.4.5分別用梯形公式和改進的歐拉公式求解區(qū)間0,2上的初值問題業(yè)=3y+8x7,y(0)=1,取精度為10,與精確解作比較，在同一個坐標(biāo)系中作出dx圖形.解輸入程序>>x0=0;y0=1;b=2;tol=1.e-1;Ht,X,

36、Yt,k,h,Pt=odtixing2(unfcn,x0,b,y0,tol)holdonS1=8/3*X-29/9+38/9*exp(-3*X),plot(X,S1,'b-')holdoffholdonH,X,Y,k,h,P=gaiEuler(unfcn,x0,b,y0,tol)holdofflegend('用梯形公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解','dy/dx=8x-3y-7,y(0)=1在0,2上的精確解','用改進的歐拉公式計算dy/dx=8x-3y-7,y(0)=1在0,2上的數(shù)值解'),

37、juwY=S1-Y;xiwY=juwY./Y;L=P,S1,juwY,xiwY運行后屏幕顯示取精度為10,，分別用梯形公式和改進的歐拉公式求此解初值問題在區(qū)間0,2上的自變量X處數(shù)值解Y,Y和精確解S及其圖形，步長HY的相對誤差xiwY和絕對誤差juwY（略）.格-庫塔方法二階龍格-庫塔方法的MATLA勰序用二階龍格一庫塔方法求解常微分方程初值問題(10.5)的數(shù)值解的MATLA駐程序functionk,X,Y,fxy,wch,wucha,P=RK2(funfcn,fun,x0,b,C,y0,h)x=x0;y=y0;p=128;n=fix(b-x0)/h);fxy=zeros(p,1);wuc

38、ha=zeros(p,1);wch=zeros(p,1);X=zeros(p,1);Y=zeros(p,length(y);k=1;X(k)=x;丫(k,:)=y'%繪圖.clc,x,h,y%十算%fxy=fxy(:);fork=2:n+1x=x+h;a2=C(3);b21=C(4);c1=C(1);c2=C(2);x1=x+a2*h;k1=feval(funfcn,x,y);y1=y+b21*h*k1;k2=feval(funfcn,x1,y1);fxy(k)=feval(fun,x);y=y+h*(c1*k1+c2*k2);X(k)=x;Y(k,:)=y;k=k+1;plot(X,

39、Y,'mh',Xfxy,'bo')grid,xlabel('自變量X'),ylabel('因變量Y')legend('用二階龍格-庫塔方法計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解，y/dx=f(x,y),y(x0)=y0的精確解y=f(x)')end%計算誤差.fork=2:n+1wucha(k)=norm(Y(k-1)-Y(k);wch(k)=norm(fxy(k)-Y(k);endX=X(1:k);Y=Y(1:k,:);fxy=fxy(1:k,:);n=1:k;wucha=wucha(

40、1:k,:);wch=wch(1:k,:);P=n',X,Y,fxy,wch,wucha;例10.5.1用二階龍格一庫塔方法求初值問題dy=1-2xy2,0_x_2,dx13xjy|x-0=0的數(shù)值解，取c1=1/4,c2=3/4,a2=b21=2/3,h-1/4,并計算與精確解的誤差,畫出精確解和數(shù)值解的圖形.解在MATLAB：作窗口輸入下面的程序>>x0=0;b=2;C=1/4,3/4,2/3,2/3;y0=0;h=1/4;k,X,Y,fxy,wch,wucha,P=RK2（unfcn,un,x0,b,C,y0,h）將運行后計算的結(jié)果列入（略），畫出精確解和數(shù)值解的圖形

41、三階龍格-庫塔方法及其MATLA醒序用三階龍格-庫塔方法求解常微分方程初值問題(10.5)的數(shù)值解的MATLA人程序functionk,X,Y,fxy,wch,wucha,P=RK3(funfcn,fun,x0,b,C,y0,h)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;丫(k,:)=y'%繪圖.clc,x,h,y%十算%fxy=fxy(:);fork=2:n+1x=x+h;a2=C(4

42、);a3=C(5);b21=C(6);b31=C(7);b32=C(8);c1=C(1);c2=C(2);c3=C(3);x1=x+a2*h;x2=x+a3*h;k1=feval(funfcn,x,y);y1=y+b21*h*k1;k2=feval(funfcn,x1,y1);y2=y+b31*h*k1+b32*h*k2;k3=feval(funfcn,x2,y2);僅y(k)=feval(fun,x);y=y+h*(c1*k1+c2*k2+c3*k3);X(k)=x;Y(k,:)=y;k=k+1;plot(X,Y,'rp',X,fxy,'bo'),grid,x

43、label('自變量X'),ylabel('因變量Y')legend('用三階龍格-庫塔方法計算dy/dx=f(x,y),y(x0)=y0在x0,b上的數(shù)值解','y/dx=f(x,y),y(x0)=y0的精確解y=f(x)')end%計算誤差.fork=2:n+1wucha(k)=norm(Y(k-1)-Y(k);wch(k)=norm(fxy(k)-Y(k);endX=X(1:k);Y=Y(1:k,:);僅丫=僅丫(1:");n=1:k;wucha=wucha(1:k,:);wch=wch(1:k,:);P=n

44、9;,X,Y,fxy,wch,wucha;例10.5.2用常用的三階龍格-庫塔公式求初值問題”一學(xué)，0*2,dx1x2yX=0的數(shù)值解，取g=1/6,c2=4/6,c3=1/6,a2=d1=1/2,a3=1,b31=1,b32=2,h=1/4,并計算與精確解的誤差，畫出精確解和數(shù)值解的圖形解在MATLAB：作窗口輸入下面的程序>>x0=0;b=2;c1=1/6;c2=4/6;c3=1/6;a2=1/2;a3=1;b21=1b31=-1;b32=2;C=c1,c2,c3,a2,a3,b21,b31,b32;y0=0;h=1/4;k,X,Y,fxy,wch,wucha,P=RK3(funfcn,fun,x0,b,C,y0,h)將運行后計算的結(jié)果(略)，畫出精確解和數(shù)值解的圖形10.5.4四階龍格-庫塔方法及其MATLA醒序用一般四階龍格-庫塔方法(10.42)式求解常微分方程初值問題(10.5)的數(shù)值解的MATLABE程序1functionk,X,Y,fxy,wch,wucha,P=RK4(funfcn,fun,x0,b,C,y0,h)x=x0;y=y0;p=128;n=fix(b-x0)/h);fxy=zeros(p,1);wucha=zeros(p,