数学11常微分方程课件

数学实验高等数学常微分方程常微分方程一、常微分方程的解析解（符号解）一、常微分方程的解析解（符号解）eq1,eqn是微分方程，其中Dy表示，Dny表示初始或边界条件 写成输出写成命令S=dsolve(eq1,eqn,独立变量)S=dsolve(eq1,eqn,初始条件,独立变量)功能求微分方程组的通解求微分方程组的特解备注独立变量的默认值为t输出S是构架数组 时，表示输出n个解。表 1 dsolve函数表例例1 求解伯努力方程并绘制积分曲线。创建dsolve01.m文件 y=dsolve(x*Dy+y=x*y2*log(x),x)%绘制积分曲线for i=-3:3 ezplot(subs(y,C1,i),0+eps,4)hold onendhold offaxis(0,4,-3,3)运行结果：y=2/x/(-log(x)2+2*C1)验算微分方程的解创建yansuan01.m文件syms x;y=dsolve(x*Dy+y=x*y2*log(x),x);L=simple(x*diff(y)+y)R=x*y2*log(x)运行结果：L=4/x*log(x)/(log(x)2-2*C1)2 R=4/x*log(x)/(-log(x)2+2*C1)2 例例2 求解柯西问题并作出积分曲线。创建dsolve01.m文件y=dsolve(y*D2y=2*(Dy2-Dy),y(0)=1,Dy(0)=2,x)ezplot(y,-3*pi/4,pi/4)运行结果：y=tan(x+1/4*pi)例例3 求解微分方程并作积分曲线。创建dsolve03.m文件syms t x y;S=dsolve(Dx+2*x+Dy+y=t,5*x+Dy+3*y=t2);%按sin(t)合并同类项S=collect(S.x,S.y,sin(t)for i=1:2 subplot(1,2,i)for j=-3:3 ezplot(subs(S(i),C1,C2,j,j+5);hold on end hold offend运行结果：S=(2*C2+3*C1)*sin(t)+C1*cos(t)+t-t2+3,(-3*C2-5*C1)*sin(t)-4+C2*cos(t)+2*t2-3*t二、常微分方程数值解二、常微分方程数值解 表表2 ode23及及ode45函数表函数表命命令令T,Y=ode23(F,tspan,Y0)T,Y=ode45(F,tspan,Y0)方方法法组合二、三阶龙格库塔算法组合二、三阶龙格库塔算法组合四、五阶龙格库塔算法组合四、五阶龙格库塔算法说说明明输输入入F是微分方程（或微分方程组）的函数文件是微分方程（或微分方程组）的函数文件tspan=t0,tn是积分区间或是积分区间或tspan=t0,t1,tn表示计算这些自表示计算这些自变量上的微分方程的解。变量上的微分方程的解。Y0是初始条件是初始条件输输出出T是自变量列向量是自变量列向量Y（一个微分方程）是因变量列向量；或（一个微分方程）是因变量列向量；或Y(微分方程组微分方程组)是矩阵，是矩阵，矩阵的行数与矩阵的行数与T的长度相等，矩阵的列数是方程的个数。的长度相等，矩阵的列数是方程的个数。例例4 用命令ode23求解初值问题并作图与解析解比较。建立函数文件fod1.mfunction F=fod1(x,y)F=-(1-2*x)/x2*y+1;创建ode01.m文件x,y=ode23(fod1,1,2,0);yy=dsolve(Dy+(1-2*x)/x2*y=1,y(1)=0,x);yy=simple(yy)err=(abs(numeric(x.2-exp(1./x).*x.2./exp(1)-y)ezplot(yy,1,2)hold onplot(x,y,r.)hold off运行结果：yy=x2-exp(1/x)*x2/exp(1)err=1.0e-004*0 0.0000 0.0000 0.0000 0.0000 0.0016 0.0353 0.0784 0.1277 0.1815 0.2389 0.2993 0.3623 0.4277 0.4954 0.5102例例 5 设一阶微分方程组当t=0时x(0)=2，y(0)=0，求微分方程在t0,50上的数值解，并画出x(t)，y(t)的解曲线以及满足微分方程的解的轨迹。步骤：（1）将微分方程组变形为 创建solve01.m文件syms x y tA=1,-1;1,1;b=-2*x+10*cos(t);-2*y+4*exp(-2*t);f=Ab 运行结果f=-x+5*cos(t)-y+2*exp(-2*t)-y+2*exp(-2*t)+x-5*cos(t)（2）建立函数文件fdo2.mfunction f=fdo2(t,y)f=-y(1)+5*cos(t)-y(2)+2*exp(-1*t);y(1)+2*exp(-2*t)-y(2)-5*cos(t);创建ode02.m文件T,Y=ode45(fdo2,0,15,2;0);figure(1)plot(T,Y(:,1),g-,T,Y(:,2),r-)legend(x(t),y(t),0)figure(2)plot(Y(:,1),Y(:,2)title(解轨迹)例例 6 求古典范得堡微分方程在初值时的数值解（t0,30）。令将微分方程转化为方程组建立函数文件vdpol.mfunction f=vdpo1(t,x)f=x(2);2*(1-x(1)2)*x(2)-x(1);创建ode02.m文件T,X=ode23(vdpo1,0,30,1;0);figure(1)plot(T,X(:,1),b-,T,X(:,2),r-)legend(x(t),dx/dt)figure(2)plot(X(:,1),X(:,2)title(函数v(x)例例7 一链条悬挂在一钉子上，起动时一端离开钉子8m，另一端离开钉子12m，摩擦为链条1m长的重量，求链条滑下来所需要的时间。设时刻t时，链条较长的一端下滑s米，并设链条的密度为，则向下拉链条的作用力由牛顿第二定理得微分方程 由此得链条下滑初值问题 用软件求解步骤（1）求初值问题符号解（位移、速度，取 计算）；（2）作位移及速度的图形，观察变化情况；（3）求链条全部滑落时的时间和速度。创建dslove04.m文件syms t s;s=dsolve(D2s=9.8/10*s-1.05*9.8,s(0)=12,Ds(0)=0)v=diff(s)subplot(1,2,1)ezplot(s,0,3)title(链条下滑位移图)subplot(1,2,2)ezplot(v,0,3)title(链条下滑速度图)f=finverse(s);t=subs(f,t,s);s=20;T=eval(t)vT=subs(v,t,T)运行结果：s=21/2+3/2*cosh(7/10*2(1/2)*t)v=21/20*sinh(7/10*2(1/2)*t)*2(1/2)T=2.5584vT=9.2865结果：（1）位移函数速度函数（2）链条滑下来需要2.5584秒，链条滑下来时的速度是9.2865m/s。上上 机机 练练 习习 题题1、用符号法求解下列微分方程，并绘制积分曲线（1）（2）；（3）（4）。2、用数值法求解下列微分方程，并作积分曲线。（1）（2）；（3）（4）3、列车以20m/s的速度行驶，当制动列车时可获得0.4m/s2的加速度。（1）求使列车停下的制动时间及制动期间列车滑行的路程。（2）作位移和速度的图形。4、由于化工生产过程中，经常要排出一些不利于环境的物质。为了保持车间内的环境卫生，必须通入大量的新鲜空气。设有一个303012m3的车间，其中的空气含有0.12%的CO2。现以含有CO2为0.04%的新鲜空气输入，问每分钟应输入多少立方米才能在10分钟后使车间空气中CO2的含量不超过0.06%?