实验六.常微分方程(组)求解.docVIP

常微分方程（组）求解 【实验目的】 1．求微分方程的符号解。 2．求微分方程的数值解。 【实验相关知识】 (一)微分方程的符号解. 微分方程的符号解也叫做解析解. 求微分方程（组）的符号解用命令dsolve. 命令格式如下： s=dsolve(方程1,方程2,…,初始条件1,初始条件2,…,自变量) 说明：用字符串表示方程，自变量缺省则默认为t. 导数用D表示, 2阶导数用D2表示,以此类推。返回值s是符号解. 例1. 求初值问题的解。 dsolve 命令格式是： dsolve(Dy = y*t/(t-5),y(0) = 2) ans = -2/3125*exp(t)*(t-5)^5 或者 s=dsolve(Dy = y*t/(t-5),y(0) = 2) s = -2/3125*exp(t)*(t-5)^5 二阶或更高阶方程的处理情况相类似。例如： , y(0) = -1, y(0) = 2 写成如下的形式输入MATLAB： dsolve(D2y - y = 0,y(0) = -1,Dy(0) = 2) ans = -3/2*exp(-t)+1/2*exp(t) 例2: 求满足初始条件的符号解. 写成如下的形式输入MATLAB： s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x) simplify(s) %如果得到符号解比较复杂,可以试试化简 pretty(s) %显示漂亮的形式。 大家得到的结果是什么呢？ (二)微分方程的数值解. 一般说来,只有对一些典型的常微分方程,才能求出它们的一般解.然而在实际问题中遇到的常微分方程往往很复杂,在许多情况下得不出一般解.所以一般是要求在若干点的近似数值解. 求数值解的命令如下: [xout,yout]=ode45(equation,[x0,xm],y0) 说明：(1) 返回值中，xout表示自变量的取值点(x0,x1,…,xn), yout表示数值解,它是一个矩阵,它的每一列对应y的一个分量。 (2) 这里equation必须是事先定义的表示微分方程(组)的M-文件。 (3) [x0,xm]是自变量的区间。 (4) y0是初始向量值。 (5)ode45还可以换成其他算法,如ode23. 注意: 命令ode45或ode23是对一阶常微分方程(组)设计的,因此对高阶常微分方程,需将它转化为一阶常微分方程组. 例如对二阶常微分方程,通过令得到一阶常微分方方程组. 在求数值解时,我们往往将数值解与画图结合, 将数值解用图像呈现出来. 例3: 求解微分方程组 先定义M-文件fun1.m。 function x=fun1(t,y) x=[y(1)*y(2)-0.3*y(1),-y(1)*y(2)] ; 然后 ts=[0,50]; % 写成ts=0:50也行 y0=[0.02,0.98]; % 表示初始值 [t,y]=ode45(fun1,ts,y0); plot(t,y(:,1),t,y(:,2)) %根据x的第一、二列同时作两条曲线 grid %为了观察方便,可添上网格线. 例4: 求解微分方程.先求符号解,再求数值解, 并作图进行比较. s=dsolve(Dy=-y+t+1,y(0)=1,t) simplify(s) 可得符号解为y=t+exp(-t). 为了求数值解,先编写M-文件fun.m function f=fun2(t,y) f=-y+t+1; 保存，再运行如下命令: clear t=0:0.02:1; y=t+exp(-t); plot(t,y) %画符号解的图形 hold on %保留已画好的图形, [t,y]=ode45(fun2,[0,1],1); plot(t,y,ro); %画数值解图形,用红色小圆圈 xlabel(t),ylabel(y) %标上各坐标名称 运行结果见下图, 可见符号解和数值解吻合得很好. 作业： 问题1：求满足初始条件的符号解。 问题 2：解微分方程组初值问题： 试画出其相图。 . 问题3：求解微分方程 先求符号解,再求数值解, 并将符号解和数值解的图形进行比较.