第9章微分方程

1第9章

MATLAB微分方程求解

29.1常微分方程初值问题常微分方程的初值问题是科学计算中的常见问题，MATLAB提供了求解常微分方程（ODE）初值问题的一系列ode函数。ode函数主要采用Runge-Kutta法求解常微分方程。常用ode函数函数名含义函数名含义ode23普通2-3阶法解ODEode23s变阶法解刚性odeode45普通4-5阶法解ODEode23t解适度刚性odeode113普通变阶法解ODEode23tb低阶法解刚性ode39.1.1ODE的机理将高阶常微分方程的初值问题，转变为1阶常微分方程组的初值问题。9.1.2普通2-3阶法解ODE例9-1编制常微分方程的函数文件dfun.mfunctiondy=dfun(x,y)dy=y-2*x/y编制如下计算程序clear;clc;X12=linspace(0,5,50);[X,Y]=ode23('dfun',X12,2);plot(X,Y)运行结果如图9-1所示图9-149.1.3普通4-5阶法解ODE例9-2编制微分方程组的函数文件dfun.mfunctiondy=dfun(x,y)dy(1,1)=-2*y(1)+y(2)+2*sin(x);dy(2,1)=10*y(1)-9*y(2)+9*(cos(x)-sin(x));编制如下程序clear;clc;X12=linspace(0,9,30);ode45('dfun',X12,[2,3]);运行结果如图9-2所示图9-259.2常微分方程边值问题9.2.1边值问题简介一般来说，微分方程边值问题可能有解、也可能无解，可能有唯一解、也可能有无数解。在假定有唯一解的前提下，边值问题有三种基本解法：迭加法：假如微分方程和边界条件是线性的，边值问题可以转化为初值问

题。试射法：将问题转化为对漏缺初值的搜索，一旦漏缺初值确定，就可以按

初值问题求解。松弛法：先猜测满足边界条件的区间网点上的解值，然后利用微分方程进

行迭代改善。69.2.2bvp4c求解思路为求解常微分方程的边值问题，MATLAB提供了bvp4c函数文件。该文件是依据有限元法中的配置法编写的。其解题的基本思路如下：

把待解问题转化为标准边值问题为期望解指定初始猜测利用原微分方程，构造残差函数用原微分方程和边界条件，借助迭代使残差不断减小为了bvp4c使用的方便，MATLAB提供了能产生较好猜测解的bvpinit指令，减轻了用户进行初始猜测的负担。此外解算指令bvp4c也能处理包含未知参数的微分方程的边值问题。79.2.3bvp4c配套指令例9-3编制微分方程组的函数文件dfun.mfunctiondy=dfun(x,y,c)dy(1,1)=y(2);dy(2,1)=-c*abs(y(1));编制残差方程组的函数文件bfun.mfunctionres=bfun(ya,yb,c)res(1,1)=ya(1);res(2,1)=yb(1)+2编制如下计算程序sol=bvp4c('dfun','bfun',sinit,[],c);plot(sol.x,sol.y(1,:),'r--')holdonc=1.5;sol=bvp4c('dfun','bfun',sinit,[],c);plot(sol.x,sol.y(1,:),'b-')legend('c=1','c=1.5')运行结果如图9-3所示clear;clc;x=linspace(0,4,5);v=[1;0];sinit=bvpinit(x,v);c=1;图9-389.3常微分方程的解析解9.3.1微分方程的通解例9-4编制如下程序clear;clc;fun='Dx=-a*t'x=dsolve(fun,'t')运行结果为fun=Dx=-a*t

x=

-1/2*a*t^2+C1例9-5编制如下程序clear;clc;fun1='Df=f+g';fun2='Dg=g-f';[f,g]=dsolve(fun1,fun2,'t')运行结果为f=

-exp(t)*(C1*cos(t)-C2*sin(t))

g=

exp(t)*(C1*sin(t)+C2*cos(t))99.3.2微分方程的特解例9-6编制如下程序clear;clc;fun='D2y=4'x=dsolve(fun,'y(0)=1','y(1)=3','x')运行结果为fun=D2y=4

x=

2*x^2+1例9-7编制如下程序clear;clc;fun1='Df=f+g';fun2='Dg=g-f';v1='f(0)=2';v2='g(0)=2';[f,g]=dsolve(fun1,fun2,v1,v2,'t')运行结果为f=

exp(t)*(2*cos(t)+2*sin(t))

g=

exp(t)*(-2*sin(t)+2*cos(t))10例9-8编制如下程序clear;clc;fun='D3y=-y';v1='y(0)=3';v2='Dy(0)=0';v3='D2y(0)=0';y=dsolve(fun,v1,v2,v3,'t')运行结果为y=

exp(-t)+2*exp(1/2*t)*cos(1/2*3^(1/2)*t)119.4偏微分方程求解9.4.1偏微分方程简介很多科学计算问题都可以归结为求解含定解条件(初始条件和边界条件)的偏微分方程问题。MATLAB提供了专门求解偏微分方程问题工具箱(PDEToolbox)。MATLAB可以求解的线性偏微分方程类型包括：

椭圆型方程抛物线方程双曲型方程特征值方程129.4.2偏微分方程的求解过程利用MATLAB求解偏微分方程的过程如下：1.完成几何造型2.输入定解条件3.选择方程类型4.划分网格5.数值求解6.结果可视化13例9-91415169.4.3PDE图形用户界面简介例9-10179.5综合实例9.5.1有阻尼自由振动例9-11编制标准一阶常微分方程组的函数文件dfun.mfunctiondz=dfun(t,z)a=-2;b=-100;dz(1,1)=z(2);dz(2,1)=a*z(2)+b*z(1);编制如下程序clear;clc;t12=linspace(0,2,100);[t,z]=ode45('dfun',t12,[1,0]);plot(t,z(:,1))gridon运行后振动规律如图9-19所示图9-19189.5.2点的运动轨迹例9-12编制标准微分方程组的函数文件dfun.mfunctiondu=dfun(t,u)du(1,1)=u(3);du(2,1)=u(4);du(3,1)=0.2*u(1);du(4,1)=0.2*u(2);编制如下程序clear;clc;t12=linspace(0,2,100);[t,u]=ode45('dfun',t12,[2,0,0,1])plot(u(:,1),u(:,2))gridon运行后，得到轨迹如图9-20所示图9-20199.5.3梁的挠曲线编制标准微分方程组的函数文件dfun.mfunctiondu=dfun(x,u)du(1,1)=u(2);du(2,1)=-0.02*x+0.01*x^3;根据定解条件,编制残差方程组的函数文件bfun.mfunctionres=bfun(ya,yb)res(1,1)=ya(1);r