Mat-lab常微分方程的求解实验报告

实验五 利用Matlab 求解常微分方程(组)的实验报告学院：数计学院 班级：1003班 姓名：黄晓丹 学号：1051020144一.实验目的:熟悉Matlab 软件中关于求解常微分方程的各种命令. 掌握利用Matlab 软件进行常微分方程的求解。

二.相关知识在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解.用字符串方程表示，自变量缺省值为t 。

导数用D 表示，2阶导数用D2表示，以此类推。

S 返回解析解。

在方程组情形，s 为一个符号结构。

[tout,yout]=ode45(‘yprime’,[t0,tf],y0) 采用变步长四阶Runge-Kutta 法和五阶Runge-Kutta-Felhberg 法求数值解，yprime 是用以表示f(t,y)的M 文件名，t0表示自变量的初始值，tf 表示自变量的终值，y0表示初始向量值。

输出向量tout 表示节点(t 0,t 1, …,t n )T ,输出矩阵yout 表示数值解，每一列对应y 的一个分量。

若无输出参数，则自动作出图形.三.实验内容:例1 求下列微分方程的解析解（1）（2）（3）方程（1）求解的MATLAB 代码为：>>clear;>>s=dsolve('Dy=a*y+b')结果为s =-b/a+exp(a*t)*C1方程（2）求解的MATLAB 代码为：>>clear;>>s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')>>simplify(s) %以最简形式显示s结果为s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x)ans =-2/3*sin(x)*cos(x)+5/3*sin(x)方程（3）求解的MATLAB 代码为：b ay y +='1)0(',0)0(,)2sin(''==-=y y y x y 1)0(',1)0(',','==-=+=g f f g g g f f>>clear;>>s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1')>>simplify(s,f) %s 是一个结构>>simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例２ 求解微分方程先求解析解，再求数值解，并进行比较。

数学与信息科学系实训报告常微分方程初值问题数值解法及其MATLAB实现姓名：学号：专业：信息与计算科学年级：2010级指导教师：朱耀生完成时间：2013年11月25—2013年12月6实验目的:(1)研究满足给定方程的可微函数的数值方法(2)培养matlab编程和上机调试能力实验基本原理和内容:根据给定的初始条件，确定常微分方程惟一解的问题叫常微分方程初值问题。

大多数实际问题难以求得解析解，必须将微分问题离散化，用数值方法求其近似解。

一阶常微分方程的初值问题的提法是，求出函数，使满足条件(1)利用数值方法解问题(1)时，通常假定解存在且惟一，解函数及右端函数具有所需的光滑程度。

数值解法的基本思想是:先取自变量一系列离散点,把微分问题(1)离散化，求出离散问题的数值解,并以此作为微分问题解的近似。

例如取步长>0,以剖分区间[,]，令=+,把微分方程离散化成一个差分方程。

以()表微分方程初值问题的解,以表差分问题的解，就是近似解的误差，称为全局误差。

因此，设计各种离散化模型,求出近似解,估计误差以及研究数值方法的稳定性和收敛性等构成了数值解法的基本内容1.欧拉法和后退欧拉法原理:function [x,y]=meuler(df,xspan,y0,h)%用途：改进欧拉公式解常微分方程y'=f(x,y), y(x0)=y0%格式：[xy]=meuler(df,a,b,y0,h) df为函数f(x,y), xspan为求解%区间[x0,xn], y0为初值y(x0), h为步长, [xy]返回节点和数值解矩阵x=xspan(1):h:xspan(2); y(1)=y0;for n=1:(length(x)-1)k1=feval(df,x(n),y(n));y(n+1)=y(n)+h*k1;k2=feval(df,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;end例题给定的初值问题y′=-y+x+2，0=<x=<1Y(0)=-1,取精确解y(x)=exp(-x)+x后退欧拉法，步长h=0.003, h=0.1求在节点k=1+0.1k (k=1,2,3……10)处的数值解若>>h=0.1;y=-1;x=1;for i=1:20k1=h*oulei_wf(x,y);k2=h*oulei_wf(x+h,y+k1);y=y+0.5*k1+0.5*k2x=x+h;z=ouleij_q(x)t=y-zendz = 4,z = 3.7000,y = -0.6150z = 1.4329, z =1.4329, t = -2.0479 z =3.7150, z =3.4435, y = -0.2571 z =1.5012, z =1.5012, t =-1.7583 z =3.4571, z = 3.2114, y =0.0763 z =1.5725, z =1.5725, t =-1.4962 z =3.2237, z = 3.0013, y =0.3876 z =1.6466, z =1.6466, t =-1.2590 z =3.0124,z =2.8112,y =0.6788z =1.7231,z =1.7231,t =-1.0444 z =2.8212,z = 2.6391,y =0.9518 z =1.8019,z =1.8019,t =-0.8501z = 2.6482,z =2.4834,y =1.2084z =1.8827,z =1.8827, t =-0.6743z =2.4916,z =2.3425, y =1.4501z =1.9653,z =1.9653, t =-0.5152实验结果的分析与研究1.对于欧拉法，步长越小，精度越高，而产生的误差越小。

实验四常微分方程的数值解法 指令：[t,y]=ode23(‘fun ’,tspan,yo) 2/3阶龙格库塔方法 [t,y]=ode45(‘fun ’,tspan,yo) 4/5阶龙格库塔方法 [t,y]=ode113(‘fun ’,tspan,yo) 高阶微分方程数值方法其中fun 是定义函数的文件名。

该函数fun 必须以为dx 输出量，以t,y 为输入量。

tspan=[t0 tfina]表示积分的起始值和终止值。

yo 是初始状态列向量。

考虑到初始条件有00d , (0)0,d d , (0)0.d SSI S S tI SI I I I tββμ⎧=-=>⎪⎪⎨⎪=-=≥⎪⎩ (5.24) 这就是Kermack 与McKendrick 的SIR 仓室模型. 方程(5.24)无法求出()S t 和()I t 的解析解.我们先做数值计算。

Matlab 代码为：function dy=rigid(t,y) dy=zeros(2,1); a=1; b=0.3;dy(1)=a*y(1).*y(2)-b*y(1); dy(2)=-a*y(1).*y(2);ts=0:.5:50; x0=[0.02,0.98];[T,Y]=ode45('rigid',ts,x0); %plot(T,Y(:,1),'-',T,Y(:,2),'*') plot(Y(:,2),Y(:,1),'b--') xlabel('s') ylabel('i')任务：1 画出i （t ），2分析各参数的影响例57：求解两点边值问题：0)5(,0)1(,32==='-''y y x y y x 。

（注意：相应的数值解法比较复杂）。

y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') ↙ y =-1/3*x^3+125/468+31/468*x^4例：用数值积分的方法求解下列微分方程 π21''2t y y -=+设初始时间t0=0；终止时间tf=3*pi ；初始条件0|',0|00====x x y y 。

实验二MATLAB数值计算常微分方程的求解引言:微分方程是描述物理问题和工程问题的重要工具，也是数学和工程学科中的重要课题。

解析方法可以求得一些简单微分方程的解析解，但对于复杂问题，常常无法得到解析解。

此时，数值求解方法成为一种有效的工具。

MATLAB是一个功能强大的数值计算软件，对于求解常微分方程组也提供了多种数值方法。

本实验将介绍MATLAB中求解常微分方程(组)的方法，并通过实例来演示这些方法的应用。

一、MATLAB的常微分方程(组)求解函数MATLAB提供了多种函数用于求解常微分方程(组)，其中最常用的函数是ode45、ode23和ode15s。

这些函数采用不同的数值方法，精度和效率也不同。

下面分别对这些函数进行简单介绍：1. ode45函数:ode45函数采用了一种自适应的步长控制算法，可以在一个时间段内自动调整步长。

这个函数的语法是：[t,y] = ode45(odefun,tspan,y0,options)其中，odefun是一个函数句柄，表示待解的常微分方程组，tspan是求解区间，y0是初始条件，options是一个结构体，包含其他的参数和选项。

2. ode23函数:ode23函数也采用了自适应的步长控制算法，但相对于ode45，它是一个简化版本，只允许使用比较简单的问题。

这个函数的语法与ode45相似：[t,y] = ode23(odefun,tspan,y0,options)3. ode15s函数:ode15s函数采用了一种隐式数值方法，适用于比较刚性（stiff）的常微分方程。

这个函数的语法与前两个函数也相似：[t,y] = ode15s(odefun,tspan,y0,options)以上是MATLAB提供的三种解常微分方程(组)的函数，根据问题的特点和求解的需求选择相应的函数。

二、实例演示在本实验中，我们将通过一个实例来演示如何使用MATLAB解常微分方程组。

假设有一个简单的线性常微分方程组：dy1/dt = -2y1 + y2dy2/dt = y1 - 2y2给定初始条件为y1(0)=1，y2(0)=0，求在t=[0,5]时，y1和y2的解。

北京理工大学珠海学院实验报告
ZHUHAI CAMPAUS OF BEIJING INSTITUTE OF TECHNOLOGY
班级2012电气2班学号xxxxxxxxx姓名陈冲指导教师张凯成绩
实验题目（实验三）常微分方程求解实验地点及时间JB501 2013/12/31（3-4节）
一、实验目的
1.掌握用程序语言来编辑函数。

2.学会用MATLAB编写RK4.m函数实现四阶Runge-kutta法。

二、实验环境
Matlab软件
三、实验内容
1、以书中第124页第11题为例编辑程序来实现计算结果。

2、使用MATLAB进行编写：
第一步：编写RK4.m函数，代码如下
第二步：利用上述函数编辑命令：（可见实验结果中的截图）
在此之前先建立一个名为ti124_11fun.m的M文件，代码如下
function z=ti124_11fun (x,y);
z=8-3*y;
再编代码：
得到结果y（0.4）=2.4654
四、实验题目
用四阶经典的龙格-库塔方法求解初值问题'83,(0)2
=-=，试取步长0.2
y y y
h=计算y的近似值，要求小数点后保留四位数字.
(0.4)
五、实验结果
五、总结
本次试验难度不大，只需注意符号别打错，实验过程无大问题。

………。

常微分方程课程实验报告实验名称 Matlab在常微分方程中的应用z =C2*exp(-2*t)+C3*exp(2*t)【3．1】作图(先做出x关于t，y关于t,z关于t的函数图像)：hold on;for c1=0:0.1:1%C1、C2和C3都大于0for c2=0:0.1:1for c3=0:0.1:1t=0:0.1:1;x=c1.*exp(2.*t)+c3.*exp(-t);y=c1.*exp(2.*t)+c2.*exp(-t)+exp(-2.*t)*c3;z=c2.*exp(2.*t)+exp(-2.*t)*c3;subplot(1,3,1),plot(t,x),legend('t~x') %x关于t的函数图像 subplot(1,3,2),plot(t,y),legend('t~y') % y关于t的函数图像 subplot(1,3,3),plot(t,z),legend('t~z') % z关于t的函数图像endendend【3.2】作图（再做出x、y、z的三维图像）：hold on;X=[];Y=[];Z=[];for c1=0:0.2:1%C1、C2和C3都大于0for c2=0:0.2:1for c3=0:0.2:1t=0:0.1:1;x=c1.*exp(2.*t)+c3.*exp(-t);y=c1.*exp(2.*t)+c2.*exp(-t)+exp(-2.*t)*c3;z=c2.*exp(2.*t)+exp(-2.*t)*c3;X=[X,x];Y=[Y,y];Z=[Z,z];endColumns 37 through 410.9753 1.0129 1.0521 1.0929 1.1353第4题：【1】编写函数M文件——cwf1.m(用于求近似解)function dy=cwf1(x,y)dy=(cos(x)-2*x*y)/(x^2-1)【2】编写脚本M文件——chang8.m[X,Y]=ode45('cwf1',[0,1],1)y1=Y’ %转置【3】则求出的近似解为（x取值为0~1）：>>y1=Columns 1 through 101.0000 0.9756 0.9524 0.9303 0.9093 0.8892 0.8701 0.8520 0.8347 0.8183Columns 11 through 200.8028 0.7880 0.7742 0.7611 0.7488 0.7374 0.7269 0.7172 0.7085 0.7008Columns 21 through 300.6941 0.6886 0.6843 0.6815 0.6802 0.6809 0.6837 0.6890 0.6976 0.7102Columns 31 through 400.7277 0.7519 0.7851 0.8319 0.8964 0.9910 1.1404 NaN NaN NaNColumn 41-Inf【4】在command窗口运行以下语句，用于求精确解：>> y=dsolve('(x^2-1)*Dy+2*x*y-cos(x)=0','y(0)=1','x') %精确解。