MATLAB作业5

MATLAB

1、 试求出下面线性微分方程的通解。

dy(t)dt55作业5

?13dy(t)dt44?64dy(t)dt33?152dy(t)dt22?176dy(t)dt?80y(t)?e?2t[sin(2t??3)?cos(3t)]&&假设上述微分方程满足已知条件，y(0)?1,y(1)?3,y(?)?2,y(0)?1,y(1)?2试

求出满足该条件的微分方程的解析解。

解： >> syms t y ;

y=dsolve(['D5y+13*D4y+64*D3y+152*D2y+176*Dy+80*y=','exp(-2*t)*(sin(2*t+pi/3)+cos(3*t))'],'y(0)=1','y(1)=3','y(pi)=2','Dy(0)=1','Dy(1)=2'); vpa(y,20) ans =

.20576131687242798354e-2*exp(-2.*t)*cos(3.*t)+.15538705805619602373e-1*exp(-2.*t)*sin(2.*t)+.76830587084294035587e-2*exp(-2.*t)*cos(2.*t)+98.159206062620455336*exp(-2.*t)*t+59.405044899367325899*exp(-2.*t)*t^3-106.24422608844727795*exp(-2.*t)*t^2-30.741892776456442810*exp(-2.*t)+.20576131687242798354e-2*exp(-2.*t)*sin(3.*t)+31.732152104579289128*exp(-5.*t)

2、 试求解下面微分方程的通解以及满足x(0)?1,x(?)?2,y(0)?0条件下的解析解。

????(t)?4x(t)?3y(t)?ex(t)?5xsin(4t) ? ?6t??2y(t)?y(t)?4x(t)?6x(t)?ecos(4t)??6t

[x,y]=dsolve('D2x+5*Dx+4*x+3*y=exp(-6*t)*sin(4*t)','2*Dy+y+4*Dx+6*x=exp(-6*t)*cos(4*t)','x(0)=1','x(pi)=2','y(0)=0')； >> vpa(x,10) ans =

0.08586166859*exp(t) - 0.057658489325149275828152894973755/(exp(7.549834435*t)^(1/4)*exp(t)^(13/4)) + (0.9469805542*exp(7.549834435*t)^(1/4))/exp(t)^(13/4) + (0.02481626654*cos(4.0*t))/exp(6.0*t) - (0.016682998530139342028763560499272*sin(4.0*t))/exp(6.0*t)

>> vpa(y,10) ans = - 0.28620556196983670815825462341309*exp(t) + 0.09045056185/(exp(7.549834435*t)^(1/4)*exp(t)^(13/4)) + (0.3018304533*exp(7.549834435*t)^(1/4))/exp(t)^(13/4) - (0.10607545320921207832043364760466*cos(4.0*t))/exp(6.0*t) + (0.0683488486*sin(4.0*t))/exp(6.0*t)

3、 试求出微分方程??y(x)?(2?1x?(x)?(1?)y1x)y(x)?xe2?5x的解析解通解，并求出满足

边界条件y(1)??,y(?)?1的解析解。

>> syms x y;

y=dsolve('D2y-(2-1/x)*Dy+(1-1/x)*y=x^2*exp(-5*x)') y =

C3*exp(t) + C2*exp((t*(x - 1))/x) + x^3/(exp(5*x)*(x - 1))

>> y=dsolve('D2y-(2-1/x)*Dy+(1-1/x)*y=x^2*exp(-5*x)','y(1)=pi','y(pi)=1') y =

(exp(t)*(exp((x - 1)/x) - x*exp((x - 1)/x) - pi*exp((pi*(x - 1))/x) - (x^3*exp((pi*(x - 1))/x))/exp(5*x) + (x^3*exp((x - 1)/x))/exp(5*x) + pi*x*exp((pi*(x - 1))/x)))/(exp(pi)*exp((x - 1)/x) - exp(1)*exp((pi*(x - 1))/x) - x*exp(pi)*exp((x - 1)/x) + x*exp(1)*exp((pi*(x - 1))/x)) - (exp((t*(x - 1))/x)*(exp(1) - x*exp(1) - pi*exp(pi) - (x^3*exp(pi))/exp(5*x) + (x^3*exp(1))/exp(5*x) + pi*x*exp(pi)))/(exp(pi)*exp((x - 1)/x) - exp(1)*exp((pi*(x - 1))/x) - x*exp(pi)*exp((x - 1)/x) + x*exp(1)*exp((pi*(x - 1))/x)) + x^3/(exp(5*x)*(x - 1))

4、 Lotka-Volterra扑食模型方程为??(t)?4x(t)?2x(t)y(t)?x?(t)?x(t)y(t)?3y(t)?y，且初值为x(0)?2,y(0)?3，

试求解该微分方程，并绘制相应的曲线。

>> syms x y t;

>> f=inline('[4*x(1)-2*x(1)*x(2); x(1)*x(2)-3*x(2)]','t','x'); >> [t,x]=ode45(f,[0,10],[2;3]);plot(t,x)

5.554.543.532.521.51012345678910

5、 是给出求解下面微分方程的MATLAB命令，

y(3)?y?ty??y?ty22?e?ty,?(0)???y(0)?2,yy(0)?0

并绘制出y(t)曲线。试问该方程存在解析解吗？选择四阶定步长Runge-Kutta算法求解该方程时，步长选择多少可以得出较好的精度，MATLAB语言给出的现成函数在速度、精度上

进行比较。

该方程的解析解不存在

>> f=inline('[x(2); x(3); -t^2*x(1)*x(2)-t^2*x(2)*x(1)^2+exp(-t*x(1))]','t','x'); [t,x]=ode45(f,[0,10],[2;0;0]); >> plot(t,x)

2.521.510.50-0.5-1-1.5-2012345678910

6、 试用解析解和数值解的方法求解下面的微分方程组

????(t)?e,x(t)??2x(t)?3x??(t)?4y?(t)?sint,y(t)?2x(t)?3y(t)?4x????5t?(0)?2x(0)?1,x?(0)?4y(0)?3,y

解析解：

>> syms t x y >>

[x,y]=dsolve('D2x=-2*x-3*Dx+exp(-5*t)','D2y=2*x-3*y-4*Dx-4*Dy-sin(t)','x(0)=1','Dx(0)=2','y(0)=3','Dy(0)=4') x =

17/(4*exp(t)) - 10/(3*exp(2*t)) + 1/(12*exp(5*t))

y =

100/(3*exp(2*t)) - 265/(16*exp(t)) - 71/(5*exp(3*t)) + 11/(48*exp(5*t)) + cos(t)/5 - sin(t)/10 + (51*t)/(4*exp(t)) 数值解：

function dx=apolloeq(t,x)

dx=[x(2);-2*x(1)-3*x(2)+exp(-5*t);x(4);2*x(1)-3*x(3)-4*x(2)-4*x(4)-sin(t)]; >> x0=[1;2;3;4];

>> [t,y]=ode45('apolloeq',[0,20],x0); >> plot(y(:,1),y(:,3))

3.532.521.510.50-0.500.20.40.60.811.21.4

7、 下面的方程在传统微分方程教程中经常被认为是刚性微分方程。使用常规微分方程解法

和刚性微分方程解法分别求解这两个微分方程的数值解，并求出解析解，用状态变量曲线比较数值求解的精度。

1??1?9y1?24y2?5cost?sint,y??3?1?y?2??24y1?51y2?9cost?sint,?3?y1(0)?y2(0)?1323(1)

解：function ydot = lorenzeq(t,y)

ydot=[9*y(1)+24*y(2)+5*cos(t)-1/3*sin(t);-24*y(1)-51*y(2)-9*cos(t)+1/3*sin(t)]; >> t_final=100; y0=[1/3;2/3];

>> [t,y]=ode45('lorenzeq',[0,t_final],y0); >> plot(t,y)

10.80.60.40.20-0.2-0.4-0.60102030405060708090100

>> opt=odeset; opt.RelTol=1e-6;

>> [t,y]=ode15s('lorenzeq',[0,t_final],y0,opt); >> plot(t,y) %刚性解法

10.80.60.40.20-0.2-0.4-0.60102030405060708090100

?1??0.1y1?49.9y2,y1(0)?1?y??2??50y2,y2(0)?2 ?y?yy3(0)?1??3?70y2?120y3,(2)解：function ydot = lorenzeq(t,y)

ydot=[-0.1*y(1)-49.9*y(2);-50*y(2); 70*y(2)-120*y(3)]; >> t_final=100; y0=[1;2;1];

[t,y]=ode45('lorenzeq',[0,t_final],y0); >> plot(t,y)

21.510.50-0.5-10102030405060708090100

>> opt=odeset; opt.RelTol=1e-6;

>> [t,y]=ode45('lorenzeq',[0,t_final],y0,opt); >> plot(t,y) %刚性解法 >> opt=odeset; opt.RelTol=1e-6;

>> [t,y]=ode45('lorenzeq',[0,t_final],y0,opt); >> plot(t,y)

21.510.50-0.5-10102030405060708090100

8、 试求出习题3中给出的微分方程边值问题数值解， 绘制出y(t)曲线，并和该习题得出

的解析解比较精度。

??u?u??0?22?x?y??9、 试用数值方法求解偏微分方程?ux?0,y?0?1,u??x?0,y?0??22y?0,x?0?0，并绘制出u函数曲面。

>> pdetool

本文来源：https://www.bwwdw.com/article/iay6.html