MATLAB_简介_7__数值法求解微分方程及微分方程组

Matlab解微分方程,除了上述的已知ODE外，还须有起始条件y0=y(x0)才能解方程式，即是在x=x0时，y(x)=y0。上述各个方程式的解析解(analyticalsolution)如下：,阮奇-库达(Runge-Kutta)方法是最通用的解ODE的方法，它可以依计算精确度的要求有低阶到高阶的各个计算式，,MATLAB应用阮奇-库达法的函数有ode23,ode45，其中ode23是同时以二阶及三阶阮奇-库达法求解，而ode45则是以四阶及五阶阮奇-库达法求解。其语法为ode23(dy,x0,xn,y0)，其中dy为ODE中的等式右边的函数（如之前介绍的），x0,xn是要解ODE的区间x0,xn的二个端点，y0是起始值(y0=y(x0)。而ode45的语法与ode23相同。,例一、要在区间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,例二、要在区间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,例三、要在区间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,例四、要在区间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计算，你可能无法察觉出其与ode23的解之间的差异，那是因为我们选的ODE代表的函数分布变化平缓，所以高阶方法就显现不出其优点。不过以二者在计算的误差上做比较，ode45的误差量级会比ode23要小。以下是几个例子：,%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的计算误差,x2,num_y2=ode45(g1,2,4,0.5);anl_y2=x2.3-7.5;%注意x2个数与x1不一定相同error_2=abs(anl_y2-num_y2)./abs(anl_y2);%ode45的计算误差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的误差的量级为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,另一个例子：%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的误差的量级为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的误差的量级为1.e-6holdoff,高阶常微分方程式,一个高阶常微分方程式可以利用变数改变(changeofvariables)方式改写成一组联立的一阶常微分方程式。例如以下的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