matlab 实验十一

实验十一微分方程（组）的解一、实验目的对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍．本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍Euler 折线法．二、实验原理1．dsolve('equ1','equ2',…)：Matlab 求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用Dy 表示y 关于自变量的一阶导数，用用D2y 表示y 关于自变量的二阶导数，依此类推．2．simplify(s)：对表达式s 使用maple 的化简规则进行化简．例如：syms xsimplify(sin(x)^2 + cos(x)^2)ans=13．[r,how]=simple(s)：由于Matlab 提供了多种化简规则，simple 命令就是对表达式s 用各种规则进行化简，然后用r 返回最简形式，how 返回形成这种形式所用的规则．例如：syms x[r,how]=simple(cos(x)^2-sin(x)^2)r = cos(2*x) how = combine4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解． 说明：(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一．(2) odefun 是显式常微分方程：⎪⎩⎪⎨⎧==00)(),(y t y y t f dt dy(3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解． (4) 要获得问题在其他指定时间点 ,210,,t t t 上的解，则令 tspan = ],,,[,210f t t t t （要求是单调的）． (5) 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver ，对于不同的ODE 问题，采用不同的Solver ．(6) 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．5．ezplot(x,y,[tmin,tmax])：符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为t 的取值范围．6．inline()：建立一个内联函数．格式：inline('expr', 'var1', 'var2',…) ，注意括号里的表达式要加引号．例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)例子：1. 几个可以直接用 Matlab 求微分方程精确解的例子： 例1：求解微分方程22x xe xy dxdy-=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2*x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4 说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为：syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x +=，解函数的图形如图 1：图1例3：求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++035y x dtdy e y x dtdx t在初始条件0|,1|00====t t y x 下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y);ezplot(x,y,[0,1.3]);axis auto微分方程的特解(式子特别长)以及解函数的图形均略．2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例4：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 的数值解，求解范围为区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y';plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000 >> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2例 5：求解描述振荡器的经典的 Ver der Pol 微分方程.7,0)0(',1)0(,0)1(222====+--μμy y y dt dy y dty d 分析：令,,121dt dx x y x ==则.)1(,1221221x x x dtdx x dt dx --==μ 先编写函数文件verderpol.m ： function xprime = verderpol(t,x) global mu;xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)]; 再编写命令文件vdp1.m ： global mu; mu = 7; y0=[1;0][t,x] = ode45('verderpol',[0,40],y0); x1=x(:,1);x2=x(:,2);plot(t,x1) 图形结果为图3．图33. 用 Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法．Euler 折线法求解的基本思想是将微分方程初值问题⎪⎩⎪⎨⎧==00)(),,(y x y y x f dxdy化成一个代数方程，即差分方程，主要步骤是用差商h x y h x y )()(-+替代微商dxdy,于是：⎪⎩⎪⎨⎧==-+)()),(,()()(00x y y x y x f h x y h x y k k k k 记)(,1k k k k x y y h x x =+=+，从而)(1h x y y k k +=+，则有1,,2,1,0).,(,),(1100-=⎪⎩⎪⎨⎧+=+==++n k y x hf y y h x x x y y k k k k k k例 6：用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=+=1)0(,22y y x y dxdy 的数值解(步长h 取0.4)，求解范围为区间[0,2]．解：本问题的差分方程为1,,2,1,0).2),( ),(,,4.0,1,021100-=⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=+====++n k y x y y x f y x hf y y h x x h y x k k k k k k （其中： 相应的Matlab 程序见附录 1．数据结果为： 0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074图形结果见图4：图4特别说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示，并和精确值比较，看看哪个更“精确”？（相应的 Matlab 程序参见附录 2）．三、实验内容1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解． >>syms x y>>y=dsolve(‘(x^2-1)*Dy+2*x*y-sin(x)=0’,’x ’) y= (-cos(x)+C1)/(x^2-1) >> pretty(y)-cos(x) + C1 ------------ 2 x - 1 2. 求微分方程x e y y y x sin 5'2''=+-的通解．>>syms x y>>y=dsolve(‘D2y-2*Dy+5*y=exp(x)*sin(x)’,’x ’)>>y=pretty(simple(y))1/3 exp(x) (3 sin(2 x) C2 + 3 cos(2 x) C1 + sin(x))3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dt dx 在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形． >> syms x y t>> [x,y]=dsolve('Dx+x+y=0','Dy+x-y=0','x(0)=1','y(0)=0','t')x =(1/2-1/4*2^(1/2))*exp(2^(1/2)*t)+(1/2+1/4*2^(1/2))*exp(-2^(1/2)*t)y =-(1/2-1/4*2^(1/2))*2^(1/2)*exp(2^(1/2)*t)+(1/2+1/4*2^(1/2))*2^(1/2)*exp(-2^(1/2)*t)-(1/2-1/4*2^(1/2))*exp(2^(1/2)*t)-(1/2+1/4*2^(1/2))*exp(-2^(1/2)*t)>> ezplot(x,y,[-1,1])>> axis auto>> 4. 分别用ode23、ode45 求上述第3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ．利用画图来比较两种求解器之间的差异．%定义函数%function dxdt=fun114ode4(t,x)dxdt=zeros(2,1);dxdt=[-x(1)-x(2);x(2)-x(1);];////////>>clc>>x0=[1,0];>> [t,x]=ode45('fun114ode4',[0,2],x0);>> subplot(2,2,1)>> plot(t,x(:,1));>> xlabel('t')>> ylabel('x')>> subplot(2,2,2);>> plot(t,x(:,2));>> xlabel('t');>> ylabel('y');>> subplot(2,2,3);>> plot3(t,x(:,1),x(:,2),'r--');>> xlabel('t0=0 tt=2');>> ylabel('x');>> zlabel('y');>> [t1,x1]=ode23('fun114ode4',[0,2],x0); >> [t,x]=ode45('fun114ode4',[0,2],x0); >> subplot(2,2,1)plot(t,x(:,1),t1,x1(:,1));xlabel('t')ylabel('x')subplot(2,2,2);plot(t,x(:,2),t1,x1(:,2));xlabel('t');ylabel('y');subplot(2,2,3);plot3(t,x(:,1),x(:,2),t1,x1(:,1),x1(:,2));xlabel('t0=0 tt=2');ylabel('x');zlabel('y');5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32y y x y y 的数值解(步长h 取0.1)，求解范围为区间[0,2]．>> h=0.1;>> y(1)=1.1;>> x=0.1:h:2;>> for n=1:19xn=x(n);yn=y(n);y(n+1)=yn+h*(yn-12*xn^2/yn^3);end>> x=[0,x];>> y=[1,y]; >> plot(x,y,'o-')6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos 'y x e y y x 的数值解(步长h 取0.1)，求解范围为区间[0,3]．四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：1,,2,1,0),()2,2()2,2(),()22(6,),(342312143211100-=⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+==++n k hL y h x f L L h y h x f L L h y h x f L y x f L L L L L h y y h x x x y y k k k k k k k k k k k k 相应的 Matlab 程序参见附录 2．试用该方法求解第5题中的初值问题． >> h=0.1;>> y(1)=1.1;>> f=sym('y-12*x^2/y^3');>> x(1)=0.1;>>for i=1:19I1=subs(f,{'x','y'},{x(i),y(i)});I2=subs(f,{'x','y'},{x(i)+h/2,y(i)+h/2*I1});I3=subs(f,{'x','y'},{x(i)+h/2,y(i)+h/2*I2});I4=subs(f,{'x','y'},{x(i)+h,y(i)+I3*h});y(i+1)=y(i)+h*(I1+2*I2+2*I3+I4)/6;end>> x=[0,x];>> y=[1,y];>> plot(x,y,'o--')7. 用ode45 方法求上述第6 题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较两者间的差异．>>edit%%函数定义function dy=fun117ode7(x,y)dy=y-12*x^2/y^3;/////y0=1;[x,y]=ode45(‘fun117ode7’,[0,2],y0);>> plot(x,y,'x-')8、分别用ode45和ode15s 求解微分方程的数值解，并进行比较 ⎩⎨⎧=-===-+1)1('',1)1(',0)1(3'4'''''223y y y x xy y x y x %%定义函数%%function dy=fun118ode8(x,y)dy=[y(2);y(3);(3*x^2+4*x*y(2)-x^2*y(3))/x^3;];//////>> x0=[0,-1,1];>> [x,y]=ode45('fun118ode8',[1,4],x0);>>[x1,y1]=ode15s('fun118ode8',[1,4],x0)；>> plot(x,y(:,1),x1,y1(:,1))五、附录附录1：(fulu1.m)clearf=sym('y+2*x/y^2');a=0;b=2;h=0.4;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))附录2：(fulu2.m)clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=0.1;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1l1=subs(f,{'x','y'},{x,y});l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))９２。