Matlab解微分方程,除了上述的已知ODE外，還須有起始條件y0=y(x0)才能解方程式，即是在x=x0時(shí)，y(x)=y0。上述各個(gè)方程式的解析解(analyticalsolution)如下：,阮奇-庫達(dá)(Runge-Kutta)方法是最通用的解ODE的方法，它可以依計(jì)算精確度的要求有低階到高階的各個(gè)計(jì)算式，,MATLAB應(yīng)用阮奇-庫達(dá)法的函數(shù)有ode23,ode45，其中ode23是同時(shí)以二階及三階阮奇-庫達(dá)法求解，而ode45則是以四階及五階阮奇-庫達(dá)法求解。其語法為ode23(dy,x0,xn,y0)，其中dy為ODE中的等式右邊的函數(shù)（如之前介紹的），x0,xn是要解ODE的區(qū)間x0,xn的二個(gè)端點(diǎn)，y0是起始值(y0=y(x0)。而ode45的語法與ode23相同。,例一、要在區(qū)間2,4解以下的ODE：%m-function,g1.mfunctiondy=g1(x,y)dy=3*x.2;x,num_y=ode23(g1,2,4,0.5);anl_y=x.3-7.5;plot(x,num_y,b:,x,anl_y,r-)title(Solutionofg1)xlabel(x),ylabel(y=f(x),grid,Y=3x2,例二、要在區(qū)間0,5解以下的ODE：%m-function,g2.mfunctiondy=g2(x,y)dy=-0.131*y;x,num_y=ode23(g2,0,5,4);anl_y=4*exp(-0.131*x);plot(x,num_y,x,anl_y,o)title(Solutionofg2)xlabel(x),ylabel(y=f(x),grid,Y=-0.13y,例三、要在區(qū)間0,2解以下的ODE：%m-function,g3.mfunctiondy=g3(x,y)dy=2*x*cos(y)2;x,num_y=ode23(g3,0,2,pi/4);anl_y=atan(x.*x+1);plot(x,num_y,b:,x,anl_y,r-)title(Solutionofg3)xlabel(x),ylabel(y=f(x),grid,Y=2xcos2y,例四、要在區(qū)間0,3解以下的ODE：%m-function,g4.mfunctiondy=g4(x,y)dy=3*y+exp(2*x);x,num_y=ode23(g4,0,3,3);anl_y=4*exp(3*x)-exp(2*x);plot(x,num_y,x,anl_y,o)title(Solutionofg4)xlabel(x),ylabel(y=f(x),gri,Y=3y+e2x,如果將上述方法改以ode45計(jì)算，你可能無法察覺出其與ode23的解之間的差異，那是因?yàn)槲覀冞x的ODE代表的函數(shù)分布變化平緩，所以高階方法就顯現(xiàn)不出其優(yōu)點(diǎn)。不過以二者在計(jì)算的誤差上做比較，ode45的誤差量級(jí)會(huì)比ode23要小。以下是幾個(gè)例子：,%m-function,g1.mfunctiondy=g1(x,y)dy=3*x.2;%m-file,odes1.m;%Solveanodeusingode23andode45clgx1,num_y1=ode23(g1,2,4,0.5);anl_y1=x1.3-7.5;error_1=abs(anl_y1-num_y1)./abs(anl_y1);%ode23的計(jì)算誤差,x2,num_y2=ode45(g1,2,4,0.5);anl_y2=x2.3-7.5;%注意x2個(gè)數(shù)與x1不一定相同error_2=abs(anl_y2-num_y2)./abs(anl_y2);%ode45的計(jì)算誤差holdonsubplot(2,2,1)plot(x1,num_y1,x1,anl_y1,o)title(ODE23solution),ylabel(y),subplot(2,2,2)plot(x1,error_y1)%注意二種方法的誤差title(ODE23error),ylabel(y)%ode23的誤差的量級(jí)為1.e-16subplot(2,2,3)plot(x2,num_y2,x2,anl_y2,o)title(ODE45solution),ylabel(y)subplot(2,2,4)plot(x1,error_y2)title(ODE45error),ylabel(y)%ode45的解沒有誤差holdoff,另一個(gè)例子：%m-function,g5.mfunctiondy=g5(x,y)dy=-y+2*cos(x);%m-file,odes1.m%Solveanodeusingode23andode45clg,x1,num_y1=ode23(g5,0,5,1);anl_y1=sin(x1)+cos(x1);error_1=abs(anl_y1-num_y1)./abs(anl_y1);x2,num_y2=ode45(g5,0,5,1);anl_y2=sin(x2)+cos(x2);error_2=abs(anl_y2-num_y2)./abs(anl_y2);holdon,subplot(2,2,1)plot(x1,num_y1,x1,anl_y1,o)title(ODE23solution),ylabel(y)subplot(2,2,2)plot(x1,error_1)%注意二種方法的誤差title(ODE23error),ylabel(y)%ode23的誤差的量級(jí)為1.e-4,subplot(2,2,3)plot(x2,num_y2,x2,anl_y2,o)title(ODE45solution),ylabel(y)subplot(2,2,4)plot(x2,error_2)title(ODE45error),ylabel(y)%ode45的誤差的量級(jí)為1.e-6holdoff,高階常微分方程式,一個(gè)高階常微分方程式可以利用變數(shù)改變(changeofvariables)方式改寫成一組聯(lián)立的一階常微分方程式。例如以下的n階方程式,以下即是解二階ODE的解法:u_prime(1)=u(1)*(1-u(2)2)-u(2);u_prime(2)=u(1);u(1)=y=y*(1-y2)-y;y=u(1);functionu_prime=eqns2(x,u)u_prime=u(1)*(1-u(2)2)-u(2);u(1);initial=00.25;x,num_y=ode23(eqns2,0,20,initial);subplot(2,1,1),plot(x,num_y(:,1)title(1stderivative