文档库 最新最全的文档下载
当前位置:文档库 › 用MATLAB解常微分方程

用MATLAB解常微分方程

用MATLAB解常微分方程
用MATLAB解常微分方程

实验四 求微分方程的解

一、问题背景与实验目的

实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解).

对微分方程(组)的解析解法(精确解),Matlab 有专门的函数可以用,本实验将作一定的介绍.

本实验将主要研究微分方程(组)的数值解法(近似解),重点介绍 Euler 折线法.

二、相关函数(命令)及简介

1.dsolve ('equ1','equ2',…):Matlab 求微分方程的解析解.equ1、equ2、…为方程(或条件).写方程(或条件)时用 Dy 表示y 关于自变量的一阶导数,用用 D2y 表示 y 关于自变量的二阶导数,依此类推.

2.simplify(s ):对表达式 s 使用 maple 的化简规则进行化简. 例如: syms x

simplify(sin(x)^2 + cos(x)^2) ans=1

3.[r,how]=simple(s):由于 Matlab 提供了多种化简规则,simple 命令就是对表达式 s 用各种规则进行化简,然后用 r 返回最简形式,how 返回形成这种形式所用的规则.

例如: syms x

[r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine

4.[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解. 说明:

(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一.

(2) odefun 是显式常微分方程:?????==0

0)()

,(y t y y t f dt dy

(3) 在积分区间 tspan =],[0f t t 上,从0t 到f t ,用初始条件0y 求解.

(4) 要获得问题在其他指定时间点Λ,210,,t t t 上的解,则令 tspan =

],,,[,210f t t t t Λ(要求是单调的).

(5) 因为没有一种算法可以有效地解决所有的 ODE 问题,为此,Matlab 提供了多种求解器 Solver ,对于不同的ODE 问题,采用不同的Solver .

(6) 要特别的是:ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序,其中:

ode23 采用龙格-库塔2 阶算法,用3 阶公式作误差估计来调节步长,具有低等的精度.

ode45 则采用龙格-库塔4 阶算法,用5 阶公式作误差估计来调节步长,具有中等的精度.

5.ezplot (x,y,[tmin,tmax]):符号函数的作图命令.x,y 为关于参数t 的符号函数,[tmin,tmax] 为 t 的取值范围.

6.inline():建立一个内联函数.格式:inline('expr', 'var1', 'var2',…) ,注意括号里的表达式要加引号.

例:Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)

三、实验内容

1. 几个可以直接用 Matlab 求微分方程精确解的例子:

例1:求解微分方程22x xe xy dx

dy

-=+,并加以验证.

求解本问题的Matlab 程序为:

syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2*x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4 说明:

(1) 行line1是用命令定义x,y 为符号变量.这里可以不写,但为确保正确性,建议写上;

(2) 行line2是用命令求出的微分方程的解:

1/2*exp(-x^2)*x^2+exp(-x^2)*C1

(3) 行line3使用所求得的解.这里是将解代入原微分方程,结果应该为0,但这里给出:

-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)

(4) 行line4 用 simplify() 函数对上式进行化简,结果为 0, 表明)(x y y =的确是微分方程的解.

例2:求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解,并画出解函数的图形.

求解本问题的 Matlab 程序为: syms x y

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')

ezplot(y)

微分方程的特解为:y=1/x*exp(x)+1/x* exp (1) (Matlab 格式),即x

e e y x

+=,

解函数的图形如图 1:

图1

例3:求微分方程组????

???=--=++035y x dt dy e y x dt dx t

在初始条件0|,1|00====t t y x 下的特解,

并画出解函数的图形.

求解本问题的 Matlab 程序为: syms x y t

[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y);

ezplot(x,y,[0,1.3]);axis auto

微分方程的特解(式子特别长)以及解函数的图形均略. 2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解).

例4:求解微分方程初值问题?????=++-=1)0(2222

y x

x y dx dy 的数值解,求解范围为区

间[0, 0.5].

fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y';

plot(x,y,'o-') >> x' ans =

0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000 >> y'

ans =

1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2.

图2

例 5:求解描述振荡器的经典的 Ver der Pol 微分方程

.7,0)0(',1)0(,0)1(22

2====+--μμy y y dt dy y dt

y d 分析:令,,121dt dx x y x =

=则.)1(,122

1221x x x dt

dx x dt dx --==μ 先编写函数文件verderpol.m :

function xprime = verderpol(t,x) global mu;

xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)]; 再编写命令文件vdp1.m : global mu; mu = 7; y0=[1;0]

[t,x] = ode45('verderpol',[0,40],y0); x1=x(:,1);x2=x(:,2); plot(t,x1)

图形结果为图3.

图3

3. 用 Euler 折线法求解

前面讲到过,能够求解的微分方程也是十分有限的.下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法.

Euler 折线法求解的基本思想是将微分方程初值问题

?????==0

0)(),

,(y x y y x f dx

dy

化成一个代数方程,即差分方程,主要步骤是用差商h

x y h x y )()(-+替代微商dx dy

,

于是:

??

?

?

?==-+)()),(,()

()(00x y y x y x f h x y h x y k k k k 记)(,1

k k k k x y y h x x =+=+,从而)(1h x y y k k +=+,则有

1,,2,1,0).

,(,

),(1

100-=???

??+=+==++n k y x hf y y h x x x y y k k k k k k Λ 例 6:用 Euler 折线法求解微分方程初值问题

??

???=+=1)0(,22

y y x y dx

dy 的数值解(步长h 取0.4),求解范围为区间[0,2].

解:本问题的差分方程为

1,,2,1,0).2),( ),(,

,4.0,1,021100-=????

????

?+=+=+====++n k y x y y x f y x hf y y h x x h y x k k k k k k Λ(其中: 相应的Matlab 程序见附录 1.

数据结果为:

0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074

图形结果见图4:

图4

特别说明:本问题可进一步利用四阶 Runge-Kutta 法求解,读者可将两个结果在一个图中显示,并和精确值比较,看看哪个更“精确”?(相应的 Matlab 程序参见附录 2).

四、自己动手

1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解.

2. 求微分方程x e y y y x sin 5'2''=+-的通解.

3. 求微分方程组

??????

?=-+=++00y x dt

dy y x dt

dx

在初始条件0|,1|00====t t y x 下的特解,并画出解函数()y f x =的图形. 4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解),求解区间为[0,2]t ∈.利用画图来比较两种求解器之间的差异.

5. 用 Euler 折线法求解微分方程初值问题

??

??

?=-=1)0(,

12'32y y x y y 的数值解(步长h 取0.1),求解范围为区间[0,2].

6. 用四阶 Runge-Kutta 法求解微分方程初值问题

?

?

?=-=1)0(,

cos 'y x e y y x 的数值解(步长h 取0.1),求解范围为区间[0,3].

四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法):

1,,2,1,0)

,()

2,2()

2,2()

,()

22(6,),(342312143211100-=??

???

????

?????

???

++=++=++==++++=+==++n k hL y h x f L L h y h x f L L h y h x f L y x f L L L L L h

y y h x x x y y k k k k k k k k k k k k Λ 相应的 Matlab 程序参见附录 2.试用该方法求解第5题中的初值问题. 7. 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解(近似解),从而利用画图来比较两者间的差异.

五、附录

附录 1:(fulu1.m)

clear

f=sym('y+2*x/y^2'); a=0; b=2; h=0.4;

n=(b-a)/h+1; x=0; y=1;

szj=[x,y];

for i=1:n-1

y=y+h*subs(f,{'x','y'},{x,y});

x=x+h;

szj=[szj;x,y];

end

szj

plot(szj(:,1),szj(:,2))

附录2:(fulu2.m)

clear

f=sym('y-exp(x)*cos(x)');

a=0;

b=3;

h=0.1;

n=(b-a)/h+1;

x=0;

y=1;

szj=[x,y];

for i=1:n-1

l1=subs(f,{'x','y'},{x,y});

l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});

l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});

l4=subs(f,{'x','y'},{x+h,y+l3*h});

y=y+h*(l1+2*l2+2*l3+l4)/6;

x=x+h;

szj=[szj;x,y];

end

szj

plot(szj(:,1),szj(:,2))

用MATLAB解常微分方程

实验四 求微分方程的解 一、问题背景与实验目的 实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解). 对微分方程(组)的解析解法(精确解),Matlab 有专门的函数可以用,本实验将作一定的介绍. 本实验将主要研究微分方程(组)的数值解法(近似解),重点介绍 Euler 折线法. 二、相关函数(命令)及简介 1.dsolve ('equ1','equ2',…):Matlab 求微分方程的解析解.equ1、equ2、…为方程(或条件).写方程(或条件)时用 Dy 表示y 关于自变量的一阶导数,用用 D2y 表示 y 关于自变量的二阶导数,依此类推. 2.simplify(s ):对表达式 s 使用 maple 的化简规则进行化简. 例如: syms x simplify(sin(x)^2 + cos(x)^2) ans=1 3.[r,how]=simple(s):由于 Matlab 提供了多种化简规则,simple 命令就是对表达式 s 用各种规则进行化简,然后用 r 返回最简形式,how 返回形成这种形式所用的规则. 例如: syms x [r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine 4.[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解. 说明: (1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一. (2) odefun 是显式常微分方程:?????==0 0)() ,(y t y y t f dt dy (3) 在积分区间 tspan =],[0f t t 上,从0t 到f t ,用初始条件0y 求解.

用Matlab软件求常微分方程的解或通解

《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 1、求微分方程的特解. 1、?? ?? ?===+-10)0(,6)0(034'2 2y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告 情况试高中

2、?? ???===++0)0(,2)0(044'2 2y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)

MATLAB求解常微分方程数值解

利用MATLAB求解常微分方程数值解

目录 1. 内容简介 (1) 2. Euler Method(欧拉法)求解 (1) 2.1. 显式Euler法和隐式Euler法 (2) 2.2. 梯形公式和改进Euler法 (3) 2.3. Euler法实用性 (4) 3. Runge-Kutta Method(龙格库塔法)求解 (5) 3.1. Runge-Kutta基本原理 (5) 3.2. MATLAB中使用Runge-Kutta法的函数 (7) 4. 使用MATLAB求解常微分方程 (7) 4.1. 使用ode45函数求解非刚性常微分方程 (8) 4.2. 刚性常微分方程 (9) 5. 总结 (9) 参考文献 (11) 附录 (12) 1. 显式Euler法数值求解 (12) 2. 改进Euler法数值求解 (12) 3. 四阶四级Runge-Kutta法数值求解 (13) 4.使用ode45求解 (14)

1.内容简介 把《高等工程数学》看了一遍,增加对数学内容的了解,对其中数值解法比较感兴趣,这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程,却对着方程无从下手,无法得出精确结果实在是让人难受的一件事情。 实际问题中更多遇到的是利用数值法求解偏微分方程问题,但考虑到先从常微分方程下手更为简单有效率,所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思,因此文章中不追求非常严格地证明,而是偏向如何利用工具实际求解出常微分方程的数值解,力求将课程上所学的知识真正地运用到实际方程的求解中去,在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。 文中求解过程中用到MATLAB进行数值求解,主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的,对各种方法进行MATLAB编程求解,并将求得的数值解与精确解对比,其中源程序在附录中。最后考察MATLAB中各个函数的适用范围,当遇到实际工程问题时能够正确地得到问题的数值解。 2.Euler Method(欧拉法)求解 Euler法求解常微分方程主要包括3种形式,即显式Euler法、隐式Euler法、梯形公式法,本节内容分别介绍这3种方法的具体内容,并在最后对3种方法精度进行对比,讨论Euler法的实用性。 本节考虑实际初值问题 使用解析法,对方程两边同乘以得到下式

Matlab求解微分方程(组)及偏微分方程(组)

第四讲 Matlab 求解微分方程(组) 理论介绍:Matlab 求解微分方程(组)命令 求解实例:Matlab 求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:解析解法和数值解法. 一.相关函数、命令及简介 1.在Matlab 中,用大写字母D 表示导数,Dy 表示y 关于自变量的一阶导数,D2y 表示y 关于自变量的二阶导数,依此类推.函数dsolve 用来解决常微分方程(组)的求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解. 注意,系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法,也称为常微分方程的符号解.但是,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰富的函数,我们将其统称为solver ,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解,则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题,为此,Matlab 提供

(完整版)实验七用matlab求解常微分方程

实验七 用matlab 求解常微分方程 一、实验目的: 1、熟悉常微分方程的求解方法,了解状态方程的概念; 2、能熟练使用dsolve 函数求常微分方程(组)的解析解; 3、能熟练应用ode45\ode15s 函数分别求常微分方程的非刚性、刚性的数值解; 4、掌握绘制相图的方法 二、预备知识: 1.微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数,称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F Λ 如果未知函数是多元函数,成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的,称为线性常微分方程,一般表示为 ) ()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--Λ 若上式中的系数 n i t a i ,,2,1),(Λ=均与t 无关,称之为常系数。 2.常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解,再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化,已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y Λ 设 ) 1(21,,',-===n n y y y y y y Λ,可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y ΛΛ 反过来,在许多情况下,一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的,一阶常系数线性微分方程组也可用特征根法求解。 3.微分方程的数值解法 除常系数线性微分方程可用特征根法求解,少数特殊方程可用初等积分法求解外,大部分微分方程无限世界,应用中主要依靠数值解法。考虑一阶常微分方程初值问题

用Matlab件求常微分方程解(或通解)

用Matlab件求常微分方程解(或通解)

————————————————————————————————作者:————————————————————————————————日期:

《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 一、求微分方程的特解. 1、?? ???===+-10)0(,6)0(034'22y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告

2、?? ???===++0)0(,2)0(044'22y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)

matlab求解常微分方程.docx

用matlab求解常微分方程 在MATLAB中,由函数dsolve()解决常微分方程(组)的求解问题,其具体格式如下:r 二dsolve('eql,eq2,???字condl,cond2,?.?;V) 匕ql,eq2,???*为微分方程或微分方程组,,condl,cond2,.??;是初始条件或边界条件,P是独立变量,默认的独立变量是讥 函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解。 dy _1 例1:求解常微分方程莎一石的MATLAB程序为: dsolve(* Dy=l/(x+y) 1r!x1), 注意,系统缺省的自变量为t,因此这里要把自变量写明。其中:Y=lambertw(X)表示函数关系Y*exp(Y)二X。 例2:求解常微分方程E'-y— 0的MATLAB程序为: Y2=dsolve(1y*D2y-Dy A2=01, 1x f) Y2=dsolve(!D2y*y-Dy A2=0 J )

我们看到有两个解,其中一个是常数0。 dx 心 ? —+ 5x + y = e dt 空_兀_3『= g2f 例3:求常微分方程组〔力 ' 通解的MATLAB 程序为: [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111 ) [X,Y]=dsolve(1 Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(- 2*t) T ,f x(0)=2,y(0)=0f ,f t T ) 以上这些都是常微分方程的精确解法,也称为常微分方程的符号解。但是,我们知 道,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析 解,此吋,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰 富的函数,我们将其统称为solver,其一般格式为: i°cosr, 7 =2 例4:求常微分方程组 y = 0 z 通解的MATLAB 程序为:

MATLABEuler法解常微分方程

Euler法解常微分方程 Euler法解常微分方程算法: Step 1 分别取积分上限、积分下限、步长 Step 2计算判断是否成立,成立转到Step 3,否则继续进行Step 4 Step 3 计算 Step 4 Euler法解常微分方程算程序: function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

Step 3 (1)做显性Euler预测 (2)将带入 Step 4计算判断是否成立,成立返回Step 3,否则继续进行Step 5 Step 5 改进Euler法解常微分方程算程序: function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' y=y' 例:用改进Euler法计算下列初值问题(取步长h=0.25) 输入:fun=inline('-x*y^2') gaijineuler2(fun,2,[0 5],0.25) 得到: x = 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500

matlab简介(解常微分方程绘制函数图像)

MATLAB简介 MATLAB是矩阵实验室(Matrix Laboratory)的简称,是美国MathWorks公司出品的商业数学软件,用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。 一、基本功能 MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案,并在很大程度上摆脱了传统非交互式程序设计语言(如C、Fortran)的编辑模式,代表了当今国际科学计算软件的先进水平。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等,主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵,它的指令表达式与数学、工程中常用的形式十分相似,故用MATLAB来解算问题要比用C,FORTRAN等语言完成相同的事情简捷得多,并且MATLAB也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C,FORTRAN,C++,JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用,此外许多的MATLAB爱好者都编写了一些经典的程序,用户可以直接进行下载就可以用。 二、特点 1) 高效的数值计算及符号计算功能,能使用户从繁杂的数学运算分析中解脱出来; 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的用户界面及接近数学表达式的自然化语言,使学者易于学习和掌握; 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。 三、优势 1.友好的工作平台编程环境 MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件,其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级,MATLAB的用户界面也越来越精致,更加接近Windows的标准界面,人机交互性更强,操作更简单。 2.强大的科学计算机数据处理能力 MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数,可以方便的实现用户所需的各种计算功能,可以用它来代替底层编程语言,如C和C++ 。在计算要求相同的情况下,使用MATLAB的编程工作量会大大减少。

用matlab求解常微分方程

实验六 用matlab 求解常微分方程 1.微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数,称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F 如果未知函数是多元函数,成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的,称为线性常微分方程,一般表示为 )()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++-- 若上式中的系数n i t a i ,,2,1),( =均与t 无关,称之为常系数。 2.常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解,再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化,已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y 设)1(21,,',-===n n y y y y y y ,可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y 反过来,在许多情况下,一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的,一阶常系数线性微分方程组也可用特征根法求解。 3.微分方程的数值解法 除常系数线性微分方程可用特征根法求解,少数特殊方程可用初等积分法求解外,大部分微分方程无限世界,应用中主要依靠数值解法。考虑一阶常微分方程初值问题 ???=<<=000)()),(,()('y t y t t t t y t f t y f

Matlab解微分方程(ODE+PDE)

常微分方程: 1 ODE解算器简介(ode**) 2 微分方程转换 3 刚性/非刚性问题(Stiff/Nonstiff) 4 隐式微分方程(IDE) 5 微分代数方程(DAE) 6 延迟微分方程(DDE) 7 边值问题(BVP) 偏微分方程(PDEs)Matlab解法 偏微分方程: 1 一般偏微分方程组(PDEs)的命令行求解 2 特殊偏微分方程(PDEs)的PDEtool求解 3 陆君安《偏微分方程的MATLAB解法 先来认识下常微分方程(ODE)初值问题解算器(solver) [T,Y,TE,YE,IE] = odesolver(odefun,tspan,y0,options) sxint = deval(sol,xint) Matlab中提供了以下解算器: 输入参数: odefun:微分方程的Matlab语言描述函数,必须是函数句柄或者字符串,必须写成Matlab

规范格式(也就是一阶显示微分方程组),这个具体在后面讲解 tspan=[t0 tf]或者[t0,t1,…tf]:微分变量的范围,两者都是根据t0和tf的值自动选择步长,只是前者返回所有计算点的微分值,而后者只返回指定的点的微分值,一定要注意对于后者tspan必须严格单调,还有就是两者数据存储时使用的内存不同(明显前者多),其它没有任何本质的区别 y0=[y(0),y’(0),y’’(0)…]:微分方程初值,依次输入所有状态变量的初值,什么是状态变量在后面有介绍 options:微分优化参数,是一个结构体,使用odeset可以设置具体参数,详细内容查看帮助 输出参数: T:时间列向量,也就是ode**计算微分方程的值的点 Y:二维数组,第i列表示第i个状态变量的值,行数与T一致 在求解ODE时,我们还会用到deval()函数,deval的作用就是通过结构体solution计算t 对应x值,和polyval之类的很相似! 参数格式如下: sol:就是上次调用ode**函数得道的结构体解 xint:需要计算的点,可以是标量或者向量,但是必须在tspan范围内 该函数的好处就是如果我想知道t=t0时的y值,不需要重新使用ode计算,而直接使用上次计算的得道solution就可以 [教程] 微分方程转换为一阶显示微分方程组方法 好,上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了,下面我们就具体开始介绍如何使用上面的知识吧! 现实总是残酷的,要得到就必须先付出,不可能所有的ODE一拿来就可以直接使用,因此,在使用ODE解算器之前,我们需要做的第一步,也是最重要的一步,借助状态变量将微分

matlabeuler法解常微分方程

Euler 法解常微分方程 Euler 法解常微分方程算法: Step 1 分别取积分上限、积分下限、步长 Step 2计算h n n +=判断b n ≤是否成立,成立转到Step 3,否则继续进行Step 4 Step 3 计算),(1n n n n y x hf y y +=+ Step 4 ),(1n n n n y x hf y y +=+ Euler 法解常微分方程算程序: function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

指导教师: 年 月 日 改进Euelr 法解常微分方程 改进Euler 法解常微分方程算法: Step 1 分别取积分上限、积分下限、步长 Step 2 取一个以h 为步长,a ,b 分别为左右端点的矩阵 Step 3 (1)做显性Euler 预测),(1n n i i y x hf y y +=+ (2)将1+i y 带入)],(),([2 h 111+++++=i i i i i i y x f y x f y y Step 4计算h n n +=判断b n ≤是否成立,成立返回Step 3,否则继续进行Step 5 Step 5 )],(),([2 h 111+++++=i i i i i i y x f y x f y y 改进Euler 法解常微分方程算程序: function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x'

Matlab求解微分方程(组)及偏微分方程(组).

第四讲 Matlab 求解微分方程(组) 理论介绍:Matlab 求解微分方程(组)命令 求解实例:Matlab 求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:解析解法和数值解法. 一.相关函数、命令及简介 1.在Matlab 中,用大写字母D 表示导数,Dy 表示y 关于自变量的一阶导数,D2y 表示y 关于自变量的二阶导数,依此类推.函数dsolve 用来解决常微分方程(组)的求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解. 注意,系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法,也称为常微分方程的符号解.但是,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰富的函数,我们将其统称为solver ,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解,则令 tspan 012[,,, ]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题,为此,Matlab 提供了多种求解器solver ,对于不同的ODE 问题,采用不同的solver.

Matlab求解常微分方程边值问题的方法

Matlab 求解常微分方程边值问题的方法:bvp4c 函数 常微分方程的边值问题,即boundary value problems ,简称BVP 问题,是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组(p 是未知参数),在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体,sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令;odefun 为描述微分方程组的函数文件;bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件;solinit 为一结构体,solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值,边界值条件分别对应a=solinit.x(l)和b=solinit.x(end); options 为bvpset 命令设定的可选函数,可采用系统默认值;p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式: 令1=y y ,21'=y y ,则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ; 将其写为函数文件twoode.m ; 同时写出边界条件函数对应文件twobc.m ; 分别使用结构solinit 和命令bvp4c 确定y-x 的关系; 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格,[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数,solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数(下述代码另存为工作目录下的twoode.m 文件) function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数(下述代码另存为工作目录下的twobc.m 文件) function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2];

常微分方程组的MATLAB求解方法

一、常微分方程组(ODEs) 简介 (1) 1. 简谐振动 (1) 2. 电路Vander Pol 方程 (1) 3. 生物种群的Volterra-Lotka 方程 (2) 4. 蝴蝶效应Lorenz 方程 (2) 二、MATLAB 数值求解ODEs的方法 (3) 1. 多变量常微分方程组的求解 (4) 2. 高阶常微分方程如何表示? (4) 3. 相图和极限环怎么绘制? (4) 个人在学习自动控制原理、现代控制理论、非线性动力学等课程时,经常遇到求解常微分方程组的问题。很多人知道MATLAB 是简便易行的一个工具,但是不会调用它自带的ode 求解器,往往还在自己编写单步欧拉法的程序,不仅求解精度差,而且程序不规范,还浪费了大量时间。以下我就工程中常见的一些非线性系统,利用MATLAB 自带的求解器,说明一下如何求解ODE 方程组、以及如何绘制相轨迹和极限环的问题。供相关专业工科大学生参考和借鉴。 一、常微分方程组(ODEs) 简介 以下列出了一些较为著名的非线性动力学系统的数学表达式,大都是由常微分方程组表达的。这种形式在工程中应用非常广泛,如力学中的非线性振动、航天领域的弹道计算、控制工程中的非线性系统等,由于自然界的大多数现象都表现出非线性,因此对于该种动力系统的研究以及微分方程的求解也具有重大的意义。以下列出一些工程应用中常见的一些由ODE 方程组所描述的动力系统。 1. 简谐振动 该式是一个2 阶非线性常微分方程。 2. 电路Vander Pol 方程

Fig 1. VanderPol 系统时域响应 3. 生物种群的 Volterra-Lotka 方程 Fig 2. Volterra-Lotka 方程时域响应(左) Fig 3 非线性动力学方程的极限环(右) 左图的捕食者 -猎物随时间变化的曲线表现出强烈的非线性,而状态变量 x 、y 的 变化却呈现出一个规则的鹅卵石状。 4. 蝴蝶效应 Lorenz 方程

最新Matlab求解微分方程(组)及偏微分方程(组)

最新Matlab 求解微分方程(组)及偏微分方程(组) 理论介绍:Matlab 求解微分方程(组)命令 求解实例:Matlab 求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:解析解法和数值解法. 一.相关函数、命令及简介 1.在Matlab 中,用大写字母D 表示导数,Dy 表示y 关于自变量的一阶导数,D2y 表示y 关于自变量的二阶导数,依此类推.函数dsolve 用来解决常微分方程(组)的求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解. 注意,系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法,也称为常微分方程的符号解.但是,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰富的函数,我们将其统称为solver ,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解,则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题,为此,Matlab 提供了多种求解器solver ,对于不同的ODE 问题,采用不同的solver.

MATLAB实验报告_常微分方程数值解

专业 序号 姓名 日期 实验3 常微分方程数值解 【实验目的】 1.掌握用MA TLAB 求微分方程初值问题数值解的方法; 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:N y(i+1)=y(i)+h*f(x(i),y(i)); end figure(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'); 【运行结果如下】:

用matlab求解常微分方程

1 .微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数,称为常微分方程。常微分方程的一般形式为 如果未知函数是多元函数,成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的,称为线性常微分方程,一般表示为 若上式中的系数均与无关,称之为常系数。 2 .常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程可化为, 两边积分可得通解为.其中为任意常数.有些常微分方程可用一些技巧,如分离变量法, 积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理, 从而他们的求解可归结为求一个特解和相应齐次微分方程的通解. 一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解,再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化,已给一个阶方程设,可将上式化为一阶方程组 反过来,在许多情况下,一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的,一阶常系数线性微分方程组也可用特征根法求解。 3.微分方程的数值解法 除常系数线性微分方程可用特征根法求解,少数特殊方程可用初等积分法求解外,大部分微分方程无限世界,应用中主要依靠数值解法。考虑一阶常微分方程初值问题 其中所谓数值解法,就是寻求在一系列离散节点上的近似值称为步长,通常取为常量。最简单的数值解法是Euler 法。 Euler 法的思路极其简单:在节点出用差商近似代替导数 这样导出计算公式(称为Euler 格式) 他能求解各种形式的微分方程。Euler 法也称折线法。 Euler 方法只有一阶精度,改进方法有二阶Runge-Kutta 法、四阶Runge-Kutta 法、五阶Runge-Kutta-Felhberg 法和先行多步法等,这些方法可用于解高阶常微分方程(组)初值问题。边值问题采用不同方法,如差分法、有限元法等。数值算法的主要缺点是它缺乏物 理理解。 4 .解微分方程的MATLAB^令

相关文档
相关文档 最新文档