《MATLAB简介7数值法求解微分方程及微分方程组》由会员分享，可在线阅读，更多相关《MATLAB简介7数值法求解微分方程及微分方程组（20页珍藏版）》请在金锄头文库上搜索。
1、Matlab解微分方程解微分方程除了上述的已知除了上述的已知ODE外，还须有起始条件外，还须有起始条件y0=y(x0)才能解方程式，即是在才能解方程式，即是在x=x0时，时，y(x)=y0。上述各个方程式。上述各个方程式的解析解的解析解(analyticalsolution)如下：如下：阮奇阮奇-库达库达(Runge-Kutta)方法是最通用的解方法是最通用的解ODE的方法，它可以依计算精确度的要求有低的方法，它可以依计算精确度的要求有低阶到高阶的各个阶到高阶的各个计算式，计算式，MATLAB应用阮奇应用阮奇-库达法的函数有库达法的函数有ode23,ode45，其中，其中ode23是同时以二阶
2、及三阶阮奇是同时以二阶及三阶阮奇-库达法求解，而库达法求解，而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;
3、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),gridY=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(
4、y=f(x),gridY=-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),gridY=2xcos2y例四、要在区间例四、要在区间0,3解以下的解以下的ODE：%m-function,g4.mfunctiondy=g4(x,y)dy=3*y+exp(
5、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),griY=3y+e2x如果将上述方法改以如果将上述方法改以ode45计算，你可能无法察计算，你可能无法察觉出其与觉出其与ode23的解之间的差异，那是因为我们的解之间的差异，那是因为我们选的选的ODE代表的函数分布变化平缓，所以高阶代表的函数分布变化平缓，所以高阶方法就显现不出其优点。不过以二者在计算的误方法就显现不出其优点。不过以二者在计算的误差上做比较
6、，差上做比较，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个数与个数
7、与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(ODE4
8、5solution),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%Solveanodeusingode23andode45clgx1,num_y1=ode23(g5,0,5,1);anl_y1=sin(x1)+cos(x1);error_1=abs(anl_y1-num_y1)./abs(anl_y1);
9、x2,num_y2=ode45(g5,0,5,1);anl_y2=sin(x2)+cos(x2);error_2=abs(anl_y2-num_y2)./abs(anl_y2);holdonsubplot(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-4subplot(2,2,3)plot(x2,num_y2,x2,a
10、nl_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(
11、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(1stderivativeofy),xlabel(x),gridsubplot(2,1,2),plot(x,num_y(:,2)title(y),xlabel(x),grid3利用利用MatlabM文件形式文件形式用四阶用四阶Runge-Kutta法法计算以上计算以上初值问题初值问题的解的解1练习练习Matlab解微分方程解微分方程2用用Matlab解微分方程组解微分方程组
