用MATLAB求解常微分方程实验报告

专业 序号 姓名 日期实验3 常微分方程数值解【实验目的】1．掌握用MATLAB 求微分方程初值问题数值解的方法；2．通过实例学习微分方程模型解决简化的实际问题；3．了解欧拉方法和龙格库塔方法的基本思想。

【实验内容】用欧拉方法和龙格库塔方法求下列微分方程初值问题的数值解，画出解的图形，对结果进行分析比较22'2,(1)(01),322;(0)1',(2)(010).(0)0(0)1x y y x x y e x y y x y x y y =+⎧≤≤=--⎨=⎩⎧=-≤≤⎨==⎩精确解或【解】:手工分析怎样求解【计算机求解】:怎样设计程序流程图变量说明能否将某算法设计成具有形式参数的函数形式【程序如下】:function f=f(x,y)f=y+2*x;clc;clear;a=0;b=1; %求解区间[x1,y_r]=ode45('f',[a b],1); %调用龙格库塔求解函数求解数值解；%% 以下利用Euler 方法求解y(1)=1;N=100;h=(b-a)/N;x=a:h:b;for i=1:Ny(i+1)=y(i)+h*f(x(i),y(i));endfigure(1)plot(x1,y_r,'r*',x,y,'b+',x,3*exp(x)-2*x-2,'k-');%数值解与真解图 title('数值解与真解图');legend('RK4','Euler','真解');xlabel('x');ylabel('y');figure(2)plot(x1,abs(y_r-(3*exp(x1)-2*x1-2)),'k-');%龙格库塔方法的误差title('龙格库塔方法的误差')xlabel('x');ylabel('Error');figure(3)plot(x,abs(y-(3*exp(x)-2*x-2)),'r-')%Euler 方法的误差title('Euler 方法的误差')xlabel('x');ylabel('Error');【运行结果如下】：【结果分析】：龙格库塔方法和Euler方法求解常微分方程都能获得比较好的数值解，相比较而言龙格库塔方法的数值解的精度远远要比Euler方法的数值解的精度高。

实验五 利用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 。