如何使用MATLAB求解微分方程(组)
matlab function求解微分方程
Matlab Function求解微分方程引言微分方程是描述自然和社会现象中的变化规律的数学工具,它在许多领域中都具有重要的应用价值。
在数值计算中,利用计算机求解微分方程成为一种常用的方法。
Matlab是一款强大的科学计算软件,它提供了丰富的函数库和工具箱,可以方便地求解各种类型的微分方程。
本文将介绍如何使用Matlab Function来求解微分方程,并通过实例说明其具体应用。
Matlab Function概述Matlab Function是Matlab中用于定义函数的关键字。
函数是一段完成特定任务的代码,可以接受输入参数并返回输出结果。
在求解微分方程中,可以通过定义一个函数来描述微分方程的数学形式,并使用Matlab内置的数值求解器来求解该微分方程。
通过封装微分方程的求解过程为一个函数,可以提高代码的复用性和可读性。
求解一阶微分方程定义微分方程函数首先需要定义微分方程的函数形式。
以一阶常微分方程dy/dx=f(x, y)为例,其中f(x, y)为已知函数。
在Matlab中,可以通过以下方式定义函数:function dy = f(x, y)dy = % 根据微分方程形式计算dy/dx的表达式end在函数中,输入参数 x 和 y 表示自变量和因变量,输出参数 dy 表示微分方程的导数值。
实际使用时,需要根据具体问题自行定义 f(x, y) 的表达式。
求解微分方程定义好微分方程函数后,可以使用Matlab内置的数值求解器来求解微分方程。
以求解某一点上的导数为例,可以使用以下代码:y0 = % 指定求解点的因变量值dydx = f(x0, y0); % 调用微分方程函数求解导数值通过以上代码,可以获得求解点上的导数值。
需要注意的是,求解点的自变量值和因变量值需要根据具体问题进行设定。
求解二阶微分方程转化为一阶微分方程组对于二阶常微分方程d2y/dx2=f(x, y, dy/dx),可以通过引入新的变量z=dy/dx,将其转化为一阶微分方程组。
使用Matlab进行微分方程求解的方法
使用Matlab进行微分方程求解的方法引言微分方程是数学中非常重要的一部分,广泛应用于物理、经济、工程等领域。
对于大部分微分方程的解析解往往难以求得,而数值解法则成为了一种常用的解决手段。
Matlab作为一种强大的科学计算软件,也提供了丰富的工具和函数用于求解微分方程,本文将介绍一些常见的使用Matlab进行微分方程求解的方法。
一、数值求解方法1. 欧拉方法欧拉方法是最简单的一种数值求解微分方程的方法,它将微分方程的微分项用差分的方式进行近似。
具体的公式为:y(n+1) = y(n) + hf(x(n), y(n))其中,y(n)表示近似解在第n个点的值,h为步长,f(x, y)为微分方程的右端项。
在Matlab中使用欧拉方法进行求解可以使用ode113函数,通过设定不同的步长,可以得到不同精度的数值解。
2. 中点法中点法是较为精确的一种数值求解微分方程的方法,它的计算公式为:k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)y(n+1) = y(n) + k2中点法通过计算两个斜率的平均值来得到下一个点的值,相较于欧拉方法,中点法能提供更精确的数值解。
3. 4阶龙格库塔法龙格库塔法是一类高阶数值求解微分方程的方法,其中4阶龙格库塔法是最常用的一种。
它的计算公式为:k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)k3 = hf(x(n) + h/2, y(n) + k2/2)k4 = hf(x(n) + h, y(n) + k3)y(n+1) = y(n) + (k1 + 2k2 + 2k3 + k4)/64阶龙格库塔法通过计算多个斜率的加权平均值来得到下一个点的值,相较于欧拉方法和中点法,它的精度更高。
二、Matlab函数和工具除了可以使用以上的数值方法进行微分方程求解之外,Matlab还提供了一些相关的函数和工具,方便用户进行微分方程的建模和求解。
如何使用MATLAB求解微分方程(组)ppt课件
差,输出参数,事件等,可缺省。 9
使用ODE?时如何编 写微分方程 ?
方式一:带额外参数,使用时需对参数进行赋值
function odefun(t,x,flag,R,L,C) %用flag说明R、L、C为变 量
xdot=zeros(2,1);
%表明其为列向量
xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);
2
Where ?
工程控制
ODE
医学生理
航空航天
金融分析
3
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
4
When ?
当对问题进行建模后,有常微分方程 需要求解时。
在生物建模中,经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
5
方法 ode23
ode45
数 ode113
当无法藉由微积分技巧求 得解析解时,这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下,常 微分方程往往只能求解出其
数值解。
在数学中,刚性方程是指一 个微分方程,其数值分析的解 只有在时间间隔很小时才会稳 定,只要时间间隔略大,其解 就会不稳定。
目前很难去精确地去定义哪 些微分方程是刚性方程,但是 大体的想法是:这个方程的解
y(1)=x(2);
y1
y2
y(2)= -t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
1000
求解程序ห้องสมุดไป่ตู้键步骤
[t,y]=ode45('odefun',[3.9 4.0],[2 8])
y
MATLAB求解微分方程(实验6)微分方程求解-文档资料
Experiments in Mathematics 微分方程求解
1
实验目的 实验内容
1、学会用Matlab求简单微分方程的解析解. 2、学会用Matlab求微分方程的数值解.
1、求简单微分方程的解析解. 2、求微分方程的数值解.
实验软件 MATLAB
2
微分方程的解析解
求微分方程(组)的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
(ltest.m)
23
1、 lorenz.m function xdot=lorenz(t,x) xdot=[-8/3,0,x(2);0,-10,10;-x(2),28,-1]*x; 2、ltest.m x0=[0 0 0.1]'; [t,x]=ode45('lorenz',[0,10],x0); plot(t,x(:,1),'-',t,x(:,2),'*',t,x(:,3),'+') pause plot3(x(:,1),x(:,2),x(:,3)),grid on
选择一组状态变量
x1 y, x2 y, , xn y(n1)
x1 x2 , x2 x3,
xn f (t, x1, x2 , , xn )
15
注意
x1 x2 , x2 x3,
xn f (t, x1, x2 , , xn )
1、建立M文件函数
function xdot = fun(t,x,y) xdot = [x2(t);x3(t);…;f(t, x1(t), x2(t),…xn(t))];
x=simple(x) % 将x简化 y=simple(y) z=simple(z)
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程(组)的解析解命令:
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
结 果:u = tan(t-c)
解 输入命令:dsolve('Du=1+u^2','t')
STEP2
STEP1
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
导弹追踪问题
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹,导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶,导弹的速度是5v0,求导弹运行的曲线方程.又乙舰行驶多远时,导弹将它击中? 解法一(解析法)
由(1),(2)消去t整理得模型:
解法二(数值解)
结 果 为:x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
2、取t0=0,tf=12,输入命令: [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图
图中,y1的图形为实线,y2的图形为“*”线,y3的图形为“+”线.
如何使用MATLAB求解微分方程学习资料
如何使用MATLAB求解微分方程学习资料MATLAB是一种强大的数值计算和科学编程平台,可以用于求解微分方程和微分方程组。
在使用MATLAB求解微分方程之前,需要掌握一些基础知识,包括MATLAB的基本语法和常用的求解微分方程的技术。
下面是一些学习资料和步骤,帮助您使用MATLAB求解微分方程。
1.学习MATLAB基本语法和操作:首先,您需要学习MATLAB的基本语法和常用操作。
您可以参考MATLAB的官方文档、教程和手册,以及MATLAB的在线资源和视频教程。
这些资源可以帮助您掌握MATLAB的基本操作,建立良好的编程习惯。
2.学习求解微分方程的方法:在使用MATLAB求解微分方程之前,您需要了解一些常用的求解微分方程的方法,例如数值方法和解析方法。
数值方法包括欧拉法、龙格-库塔法、四阶龙格-库塔法等;解析方法包括分离变量法、线性微分方程的常系数齐次法和非齐次法等。
您可以参考微积分的教科书、在线资源和视频教程,掌握这些方法。
3. 使用MATLAB求解一阶微分方程:一阶微分方程是最简单的微分方程形式。
您可以首先尝试使用MATLAB求解一阶微分方程。
MATLAB提供了几个函数来求解一阶微分方程,例如ode45、ode23、ode113等。
您可以使用这些函数来解决特定的一阶微分方程,并观察结果。
可以使用plot 函数绘制微分方程的解,以获得更直观的理解。
4.使用MATLAB求解高阶微分方程:一旦您熟悉了使用MATLAB求解一阶微分方程的方法,您可以尝试使用同样的方法来求解高阶微分方程。
在求解高阶微分方程时,您需要将其转化为一组一阶微分方程。
例如,对于二阶线性微分方程,您可以引入一个新的变量来表示未知函数的导数,然后将其转化为一组一阶微分方程。
然后,您可以使用相同的求解函数来求解这组一阶微分方程。
5. 使用MATLAB求解微分方程组:对于多元微分方程组,MATLAB提供了更多的函数来求解。
例如,ode45s可以用于求解刚体动力学方程,ode23t可以用于求解刚体动力学方程。
matlab 求微分方程组数值解
matlab 求微分方程组数值解使用Matlab求解微分方程组是一种常见的数值方法。
微分方程组是描述自然界中许多现象的数学模型,它们可以用一组关于未知函数及其导数的方程来表示。
通过求解微分方程组,我们可以得到未知函数在给定条件下的数值解。
在Matlab中,求解微分方程组可以使用ode45函数。
该函数是一个常用的求解常微分方程初值问题的函数,它使用四阶龙格-库塔法(RK4)进行数值求解。
使用ode45函数求解微分方程组的步骤如下:定义微分方程组。
在Matlab中,可以使用匿名函数或函数句柄的方式定义微分方程组。
例如,对于一个二阶微分方程组:dy1/dt = f1(t, y1, y2)dy2/dt = f2(t, y1, y2)可以定义一个匿名函数:f = @(t, y) [f1(t, y(1), y(2)); f2(t, y(1), y(2))]其中,t是自变量,y是未知函数的向量。
接下来,指定求解的时间区间和初值条件。
时间区间可以通过指定起始时间和结束时间来确定。
初值条件是指在起始时间处未知函数的值。
初值条件可以通过一个向量来表示。
例如,对于一个二阶微分方程组,初值条件可以表示为一个长度为2的向量。
然后,调用ode45函数进行求解。
ode45函数的输入参数包括定义的微分方程组、时间区间和初值条件。
该函数会返回数值解和对应的时间点。
可以通过绘制图形或打印数值解来展示结果。
Matlab提供了丰富的绘图函数,可以方便地将数值解可视化。
需要注意的是,求解微分方程组时,应选择合适的数值方法和步长,以保证数值解的精度和稳定性。
对于复杂的微分方程组,可能需要进行参数调整和迭代求解,以得到满意的结果。
使用Matlab求解微分方程组是一种便捷而有效的数值方法。
通过定义微分方程组、指定时间区间和初值条件,调用ode45函数进行求解,可以得到微分方程组的数值解。
这种方法在科学研究和工程实践中具有广泛的应用,可以帮助我们更好地理解和分析自然界中的现象。
matlab用四阶龙格库塔函数求解微分方程组
一、介绍Matlab作为一种强大的科学计算软件,提供了众多函数和工具来解决微分方程组。
其中,四阶龙格库塔函数是一种常用的数值方法,用于求解常微分方程组。
本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组,并对该方法的原理和实现进行详细说明。
二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法,用于求解常微分方程组。
它是一种显式的Runge-Kutta方法,通过逐步逼近微分方程的解,在每一步使用多个中间值来计算下一步的解。
该方法通过四个中间值来计算下一步的状态,并且具有较高的精度和稳定性。
三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中,可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。
ode45函数是Matlab提供的用于求解常微分方程组的函数,可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。
1. 定义微分方程组我们需要定义要求解的微分方程组。
可以使用Matlab中的匿名函数来定义微分方程组,例如:```matlabf = (t, y) [y(2); -sin(y(1))];```其中,f是一个匿名函数,用于表示微分方程组。
在这个例子中,微分方程组是y' = y2, y2' = -sin(y1)。
2. 指定初值条件和求解区间接下来,我们需要指定微分方程组的初值条件和求解区间。
初值条件可以通过指定一个初始时刻的状态向量来完成,例如:```matlabtspan = [0, 10];y0 = [0, 1];```其中,tspan表示求解区间,y0表示初值条件。
3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。
具体的调用方式如下:```matlab[t, y] = ode45(f, tspan, y0);```其中,t和y分别表示求解的时间点和对应的状态值。
四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。
matlab求解微分方程组
matlab求解微分方程组
Matlab 是一种非常强大的工具,可以用来求解各种微分方程组。
它可以解决复杂的微分方程组,有助于我们快速获得精确的解决方案。
Matlab 提供了一系列函数,用于求解微分方程组。
其中最常用的
函数是 ode45、ode15s 和 ode23s。
它们可以用来求解常微分方程组,也可以用来求解非线性方程组。
首先,我们需要准备好微分方程组的初始条件。
然后,我们可以
使用 Matlab 的 ode45 函数来求解微分方程组。
ode45 函数可以求
解常微分方程组,它使用 Runge-Kutta 方法来求解方程组。
使用 ode45 函数求解微分方程组的步骤如下:
1. 首先,我们需要准备好微分方程组的初始条件,并将其输入到Matlab 中。
2. 然后,我们需要定义一个 Matlab 函数,用于定义微分方程组。
3. 接下来,我们可以使用 ode45 函数来求解微分方程组。
ode45 函数的第一个参数是 Matlab 函数,用于定义微分方程组;第二个参数是初始条件;第三个参数是微分方程组的解的范围。
4. 最后,我们可以使用 Matlab 的 plot 函数来绘制微分方程组的解。
Matlab 提供了很多有用的函数,可以用来求解微分方程组。
它的运算速度快,可以让我们获得更准确的解决方案。
使用 Matlab 可以节省大量的时间,提高工作效率。
Matlab求解微分方程(组)及偏微分方程(组)
第四讲 Matlab 求解微分方程(组)理论介绍:Matlab 求解微分方程(组)命令 求解实例:Matlab 求解微分方程(组)实例实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:解析解法和数值解法. 一.相关函数、命令及简介1.在Matlab 中,用大写字母D 表示导数,Dy 表示y 关于自变量的一阶导数,D2y 表示y 关于自变量的二阶导数,依此类推.函数dsolve 用来解决常微分方程(组)的求解问题,调用格式为:X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解.注意,系统缺省的自变量为t2.函数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.表1 Matlab中文本文件读写函数说明:ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶微分方程(组)的初值问题的解的Matlab常用程序,其中:ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度.ode45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,具有中等的精度.3.在matlab命令窗口、程序或函数中创建局部函数时,可用内联函数inline,inline函数形式相当于编写M函数文件,但不需编写M-文件就可以描述出某种数学关系.调用inline函数,只能由一个matlab表达式组成,并且只能返回一个变量,不允许[u,v]这种向量形式.因而,任何要求逻辑运算或乘法运算以求得最终结果的场合,都不能应用inline函数,inline函数的一般形式为:FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如:(求解F(x)=x^2*cos(a*x)-b ,a,b是标量;x是向量)在命令窗口输入:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’); g= Fofx([pi/3 pi/3.5],4,1) 系统输出为:g=-1.5483 -1.7259注意:由于使用内联对象函数inline 不需要另外建立m 文件,所有使用比较方便,另外在使用ode45函数的时候,定义函数往往需要编辑一个m 文件来单独定义,这样不便于管理文件,这里可以使用inline 来定义函数. 二.实例介绍1.几个可以直接用Matlab 求微分方程精确解的实例 例1 求解微分方程2'2x y xy xe -+=程序:syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’)例 2 求微分方程'0x xy y e +-=在初始条件(1)2y e =下的特解并画出解函数的图形.程序:syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 求解微分方程组530tdx x y e dtdy x y dt⎧++=⎪⎪⎨⎪--=⎪⎩在初始条件00|1,|0t t x y ====下的特解并画出解函数的图形.程序: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 auto2.用ode23、ode45等求解非刚性标准形式的一阶微分方程(组)的初值问题的数值解(近似解)例 4 求解微分方程初值问题2222(0)1dy y x xdx y ⎧=-++⎪⎨⎪=⎩的数值解,求解范围为区间[0,0.5].程序:fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-')例 5 求解微分方程22'2(1)0,(0)1,(0)0d y dyy y y y dt dtμ--+===的解,并画出解的图形.分析:这是一个二阶非线性方程,我们可以通过变换,将二阶方程化为一阶方程组求解.令12,,7dyx y x dtμ===,则 121221212,(0)17(1),(0)0dx x x dtdx x x x x dt⎧==⎪⎪⎨⎪=--=⎪⎩ 编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 命令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)练习与思考:M-文件vdp.m 改写成inline 函数程序? 3.用Euler 折线法求解Euler 折线法求解的基本思想是将微分方程初值问题00(,)()dyf x y dxy x y ⎧=⎪⎨⎪=⎩ 化成一个代数(差分)方程,主要步骤是用差商()()y x h y x h +-替代微商dydx,于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例 6 用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解(步长h 取0.4),求解范围为区间[0,2].分析:本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序:>> 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-1y=y+h*subs(f,{'x','y'},{x,y});%subs ,替换函数 x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))说明:替换函数subs 例如:输入subs(a+b,a,4) 意思就是把a 用4替换掉,返回 4+b ,也可以替换多个变量,例如:subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha 替换a 和2替换b ,返回 cos(alpha)+sin(2)特别说明:本问题可进一步利用四阶Runge-Kutta 法求解,Euler 折线法实际上就是一阶Runge-Kutta 法,Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩相应的Matlab 程序为:>> 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-1l1=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))练习与思考:(1)ode45求解问题并比较差异. (2)利用Matlab 求微分方程(4)(3)''20y y y -+=的解.(3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解. (4)利用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解. 提醒:尽可能多的考虑解法 三.微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准形式的一阶显式微分方程(组)问题,因此在使用ODE 解算器之前,我们需要做的第一步,也是最重要的一步就是借助状态变量将微分方程(组)化成Matlab 可接受的标准形式.当然,如果ODEs 由一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍如何将它变换成一个一阶显式微分方程组.Step 1 将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高排列.形式为:()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x x y y y y ----⎧=⎨=⎩Step 2 为每一阶微分式选择状态变量,最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y--++++========注意:ODEs 中所有是因变量的最高阶次之和就是需要的状态变量的个数,最高阶的微分式不需要给它状态变量.Step 3 根据选用的状态变量,写出所有状态变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m nm n x x x x x x x f t x x x x xx xg t x x x x +++++======练习与思考:(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩ 提示:使用符号计算函数solve 求'''',x y ,然后利用求解微分方程的方法 四.偏微分方程解法Matlab 提供了两种方法解决PDE 问题,一是使用pdepe 函数,它可以求解一般的PDEs,具有较大的通用性,但只支持命令形式调用;二是使用PDE 工具箱,可以求解特殊PDE 问题,PDEtoll 有较大的局限性,比如只能求解二阶PDE 问题,并且不能解决片微分方程组,但是它提供了GUI 界面,从复杂的编程中解脱出来,同时还可以通过File —>Save As 直接生成M 代码.1.一般偏微分方程(组)的求解(1)Matlab 提供的pdepe 函数,可以直接求解一般偏微分方程(组),它的调用格式为:sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描述函数,它必须换成标准形式:(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ 这样,PDE 就可以编写入口函数:[c,f,s]=pdefun(x,t,u,du),m,x,t 对应于式中相关参数,du 是u 的一阶导数,由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界条件描述函数,它必须化为形式:(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂==∂ 于是边值条件可以编写函数描述为:[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下边界,b 表示上边界.@pdeic 是PDE 的初值条件,必须化为形式:00(,)u x t u =,故可以使用函数描述为:u0=pdeic(x)sol 是一个三维数组,sol(:,:,i)表示i u 的解,换句话说,k u 对应x(i)和t(j)时的解为sol(i,j,k),通过sol ,我们可以使用pdeval 函数直接计算某个点的函数值.(2)实例说明 求解偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中, 5.7311.46()x x F x e e -=-且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0uu t u t t x∂===∂ 解:(1)对照给出的偏微分方程和pdepe 函数求解的标准形式,原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ %目标PDE 函数function [c,f,s]=pdefun(x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)) end(2)边界条件改写为:下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦%边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) pa=[0;ua(2)]; qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1]; end(3)初值条件改写为:1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数 function u0=pdeic(x) u0=[1;0]; end(4)编写主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); subplot(2,1,1) surf(x,t,sol(:,:,1)) subplot(2,1,2) surf(x,t,sol(:,:,2))练习与思考: This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x xπ∂∂∂=∂∂∂ This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t uu t e t xπ-∂=+=∂ 2.PDEtool 求解偏微分方程(1)PDEtool (GUI )求解偏微分方程的一般步骤在Matlab 命令窗口输入pdetool ,回车,PDE 工具箱的图形用户界面(GUI)系统就启动了.从定义一个偏微分方程问题到完成解偏微分方程的定解,整个过程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω,通过公式把Matlab 系统提供的实体模型:矩形、圆、椭圆和多边形,组合起来,生成需要的平面区域.Step 2 “Boundary 模式”定义边界,声明不同边界段的边界条件.Step 3 “PDE 模式”定义偏微分方程,确定方程类型和方程系数c,a,f,d ,根据具体情况,还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω,可以控制自动生成网格的参数,对生成的网格进行多次细化,使网格分割更细更合理.Step 5 “Solve 模式”解偏微分方程,对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,设置初始边界条件后可以求出给定时刻t 的解;对于特征值问题,可以求出给定区间上的特征值.求解完成后,可以返回到Step 4,对网格进一步细化,进行再次求解.Step 6 “View 模式”计算结果的可视化,可以通过设置系统提供的对话框,显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示.(2)实例说明用法求解一个正方形区域上的特征值问题:12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩ 正方形区域为:11,1 1.x x -≤≤-≤≤(1)使用PDE 工具箱打开GUI 求解方程(2)进入Draw 模式,绘制一个矩形,然后双击矩形,在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框(3)进入Boundary 模式,边界条件采用Dirichlet 条件的默认值(4)进入PDE 模式,单击工具栏PDE 按钮,在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后关闭对话框(5)单击工具栏的 按钮,对正方形区域进行初始网格剖分,然后再对网格进一步细化剖分一次(6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特征值区域为[-20,20](7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color、Height(3-D plot)和show mesh项,然后单击Done确认(8)单击工具栏的“=”按钮,开始求解。
matlab解微分方程组
matlab解微分方程组
MATLAB是一种强大的计算工具,能够以高效的方式处理复杂的数学问题。
由于其灵活的编程接口和拥有大量可用的函数,MATLAB可以被用于解决各种不同类型的微分方程组。
本文将介绍如何使用MATLAB 解微分方程组。
MATLAB可以利用拟牛顿发展算法,利用函数ode45来解决常微分方程组(Ordinary Differential Equations,简称ODEs)。
生成积分函数,与函数ode45耦合在一起,可以用ode45函数解ODE。
第一步,将微分方程组写成一阶形式,即:dy/dx=f(x,y),其中y为未知变量,x为变量,f(x,y) 为表达式。
第二步,使用MATLAB编程生成函数解微分方程组。
函数ode45是MATLAB中用于解ODE的函数,它使用拟牛顿发展算法,可以得到非线性ODE的数值解。
首先写出解ODE的函数,接受自变量x和因变量y 做参数,并返回相应的函数值;然后,可以调用函数ode45来解这些ODE,函数将接受积分端点、积分步长和积分函数作为参数,并返回结果。
最后,将结果可视化展示出来。
使用数据可视化函数,如plot,可以将结果以曲线的形式展示出来,方便对结果进行后续处理。
总结起来,通过使用MATLAB的ode45函数,配合编写的解ODE 函数,可以快捷高效地解决一般微分方程组问题。
通过可视化函数,还可以将解决出的结果展示出来,为数据分析提供便利。
matlab解二阶微分方程组
matlab解二阶微分方程组Matlab是最为常用的数学计算软件之一,它不仅可以解一元或多元方程,还可以解微分方程。
其中,解二阶微分方程组是比较常见的问题,因此本文将对Matlab如何解二阶微分方程组进行详细介绍。
一、Matlab解二阶微分方程组的基本方法解二阶微分方程组需要用到ode45(ode23、ode113等)函数,该函数是Matlab求解微分方程组的函数之一。
它实现了一种常见的数值积分算法,用于求解Ordinary Differential Equations (ODEs)。
解二阶微分方程组的一般形式为:y''(t) = f(t,y(t),y'(t))y(a) = y0, y'(a)=y1其中,y''(t)表示时间t时的二阶导数,y(t)表示未知函数,f(t,y(t),y'(t))为已知函数,a为初值点,y0和y1为该点处的函数值和函数的一阶导数值。
二、使用Matlab解二阶微分方程组的步骤1. 用matlab的Function Handle定义所求解的方程组可使用函数句柄将所求解的方程组转化为Matlab可以理解的形式。
任何匿名函数都可以使用函数句柄定义,这样就可以将其输入到ODE求解器中。
例如,假设要解决以下方程组:y''1(t) = -2*y'1(t) + y2(t)y''2(t) = -y1(t) - 2*y'2(t)y(0) = [1 0]y'(0) = [0 1]定义该方程组的函数句柄如下:f = @(t,y) [y(2); -y(1)-2*y(2); y(4); -2*y(3)+y(4)]2. 调用ode函数进行求解ode函数有很多种,如ode45,ode23,ode113等。
其基本用法如下:[T,Y] = ode45(fun,tspan,y0);其中,fun是一个函数句柄,指向所求的函数;tspan是t的取值区间,y0是初值条件,T和Y是求解得到的值。
matlab dsolve 微分方程组
在MATLAB中,可以使用`dsolve`函数来求解微分方程组。
`dsolve`函数可以求解常微分方程(Ordinary Differential Equations,ODE)和偏微分方程(Partial Differential Equations,PDE)。
下面是一个示例,演示如何使用`dsolve`函数来求解一个简单的微分方程组:
```matlab
syms t x(t) y(t)
eq1 = @(t,x) x(t)/x(t-1) - 2; 第一个方程
eq2 = @(t,x) x(t-1)/x(t) - 3; 第二个方程
sol = dsolve({eq1, eq2}, x(t), t); 求解微分方程组
disp(sol); 显示解
```
在这个示例中,我们定义了两个方程`eq1`和`eq2`,然后使用`dsolve`函数来求解这两个方程组成的微分方程组。
注意,我们需要将方程以函数的形式传递给`dsolve`函数。
在`dsolve`函数中,第一个参数是一个包含所有方程的向量,第二个参数是要求解的未知函数。
`dsolve`函数将返回一个包含所有解的表达式。
在本例中,我们将解存储在`sol`变量中,并使用`disp`函数显示解。
请注意,在使用`dsolve`函数时,需要确保输入的方程是正确的,并且与所求解的问题相对应。
此外,还需要注意符号和函数的定义和使用方式。
matlab解微分方程组拉式变换后的方程
一、概述微分方程组是描述自然界中众多物理现象的重要数学工具,它在工程、物理学、生物学等领域有着广泛的应用。
Matlab作为一种高效的科学计算软件,能够方便地对微分方程组进行求解和分析。
在求解微分方程组时,拉普拉斯变换是一种非常重要的方法,它可以将微分方程转化为代数方程,从而简化求解的过程。
本文将重点探讨利用Matlab对微分方程组进行拉普拉斯变换后的求解过程。
二、微分方程组的拉普拉斯变换1. 拉普拉斯变换的定义拉普拉斯变换是一种用来处理微分方程的常用方法。
对于一个函数f(t),它的拉普拉斯变换定义如下:L{f(t)} = F(s) = ∫[0,∞]e^(-st)f(t)dt其中,s为复变量,t为实变量。
通过拉普拉斯变换,可以将微分方程组转化为代数方程组,从而利用代数方法进行求解。
2. 微分方程组的拉普拉斯变换考虑一个n阶线性微分方程组:a_n(t)y^n + a_(n-1)(t)y^(n-1) + ... + a_1(t)y' + a_0(t)y = f(t)通过拉普拉斯变换,将上述微分方程组转化为代数方程组:A_n(s)Y(s) + A_(n-1)(s)Y(s) + ... + A_1(s)Y(s) + A_0(s)Y(s) = F(s)其中,Y(s)和F(s)分别为y(t)和f(t)的拉普拉斯变换,A_n(s)等为对应的系数的拉普拉斯变换。
三、Matlab求解拉普拉斯变换后的微分方程组1. 预处理在利用Matlab求解微分方程组之前,需要对微分方程组进行拉普拉斯变换。
假设有一个简单的n阶线性微分方程组:a_n(t)y^n + a_(n-1)(t)y^(n-1) + ... + a_1(t)y' + a_0(t)y = f(t)通过拉普拉斯变换,可以得到:A_n(s)Y(s) + A_(n-1)(s)Y(s) + ... + A_1(s)Y(s) + A_0(s)Y(s) = F(s)将上述代数方程组输入Matlab中进行求解。
matlab求解二元二阶微分方程组
matlab求解二元二阶微分方程组求解二元二阶微分方程组是数学和工程领域中的重要问题,而使用Matlab可以方便地进行求解。
本文将介绍如何使用Matlab求解二元二阶微分方程组,并通过实例说明具体的求解过程。
让我们来回顾一下二元二阶微分方程组的一般形式。
设有二元函数x(t)和y(t),则二元二阶微分方程组可以写为:x''(t) = f(t, x(t), y(t))y''(t) = g(t, x(t), y(t))其中f(t, x(t), y(t))和g(t, x(t), y(t))是已知的函数。
我们的目标是求解x(t)和y(t)。
在Matlab中,可以使用ode45函数来求解二元二阶微分方程组。
该函数是Matlab中常用的求解常微分方程的函数,可以自动选择合适的步长进行求解,具有较高的准确性和稳定性。
接下来,我们将通过一个具体的例子来演示如何使用Matlab求解二元二阶微分方程组。
假设我们要求解以下二元二阶微分方程组:x''(t) = -x(t) + 2y(t)y''(t) = 3x(t) - y(t)我们需要将该方程组转化为一阶方程组。
引入新的变量u1 = x(t),u2 = x'(t),u3 = y(t),u4 = y'(t),则原方程组可以写为:u1' = u2u2' = -u1 + 2u3u3' = u4u4' = 3u1 - u3接下来,我们可以使用Matlab编写代码来求解该方程组。
首先,定义一个匿名函数,表示方程组的右端项:f = @(t,u) [u(2); -u(1) + 2*u(3); u(4); 3*u(1) - u(3)];然后,设置求解的时间范围和初始条件,并调用ode45函数进行求解:tspan = [0 10]; % 求解的时间范围u0 = [1; 0; 2; 0]; % 初始条件[t, u] = ode45(f, tspan, u0);我们可以通过绘图来显示求解结果。
matlab求解二元微分方程组
matlab求解二元微分方程组使用Matlab求解二元微分方程组是一种常见的数学问题求解方法。
二元微分方程组是指含有两个未知函数的微分方程组,通常形式为:dx/dt = f(x, y)dy/dt = g(x, y)其中,x和y是未知函数,t是自变量,f和g是给定的函数。
求解二元微分方程组的目标是找到满足上述方程组的函数x(t)和y(t)。
Matlab是一种功能强大的数值计算软件,它提供了各种数学工具和函数库,可以方便地求解二元微分方程组。
下面将介绍使用Matlab 求解二元微分方程组的步骤和方法。
我们需要定义给定的函数f和g,并确定初始条件。
在Matlab中,可以使用匿名函数来定义f和g。
例如,假设我们要求解的二元微分方程组为:dx/dt = -x + ydy/dt = -2y我们可以定义f和g如下:f = @(x, y) -x + y;g = @(x, y) -2*y;接下来,我们需要确定初始条件。
初始条件是指在某个特定的时间点上,未知函数x和y的值。
在Matlab中,可以使用初值问题来求解二元微分方程组。
初值问题需要给定初始时间点t0和初始条件x0和y0。
假设我们要求解的二元微分方程组的初始条件为:t0 = 0;x0 = 1;y0 = 2;我们可以使用ode45函数来求解二元微分方程组。
ode45函数是Matlab中常用的求解常微分方程的函数,它使用了四阶龙格-库塔方法。
使用ode45函数求解二元微分方程组的语法如下:[t, sol] = ode45(@(t, y) [f(y(1), y(2)); g(y(1), y(2))], [t0, tf], [x0; y0]);其中,t是时间点的向量,sol是解向量。
@(t, y) [f(y(1), y(2)); g(y(1), y(2))]表示一个匿名函数,它将t和y作为输入,返回[f(y(1), y(2)); g(y(1), y(2))]。
[t0, tf]是时间范围,[x0; y0]是初始条件。
MATLAB微分方程几种求解方法及程序
第五章 控制系统仿真§5.2 微分方程求解方法以一个自由振动系统实例为例进行讨论。
如下图1所示弹簧-阻尼系统,参数如下: M=5 kg, b=1 N.s/m, k=2 N/m, F=1NF图1 弹簧-阻尼系统假设初始条件为:00=t 时,将m 拉向右方,忽略小车的摩擦阻力,m x 0)0(= s m x /0)0(=•求系统的响应。
)用常微分方程的数值求解函数求解包括ode45、ode23、ode113、ode15s 、ode23s 等。
wffc1.m myfun1.m一、常微分方程的数值求解函数ode45求解 解:系统方程为 F kx x b x m =++•••这是一个单变量二阶常微分方程。
将上式写成一个一阶方程组的形式,这是函数ode45调用规定的格式。
令: x x =)1( (位移))1()2(••==x x x (速度) 上式可表示成:⎥⎦⎤⎢⎣⎡--=⎥⎦⎤⎢⎣⎡=⎥⎥⎦⎤⎢⎢⎣⎡••)1(*20)2(*101)2()2()2()1(x x x x x x x 下面就可以进行程序的编制。
%写出函数文件myfun1.mfunction xdot=myfun1(t,x)xdot=[x(2);1-10*x(2)-20*x(1)];% 主程序wffc1.mt=[0 30];x0=[0;0];[tt,xx]=ode45(@myfun1,t,x0);plot(tt,yy(:,1),':b',tt,yy(:,2),'-r') legend('位移','速度')title('微分方程的解 x(t)')二、方法2:F kx x b x m =++•••251)()()(2++==s s s F s X s G%用传递函数编程求解ksys1.mnum=1;den=[5 1 2];%printsys(num,den)%t=0:0.1:10;sys=tf(num,den);figure(1)step(sys)figure(2)impulse(sys)figure(3)t=[0:0.1:10]';ramp=t;lsim(sys,ramp,t);figure(4)tt=size(t);noise=rand(tt,1);lsim(sys,noise,t)figure(5)yy=0.1*t.^2;lsim(num,den,yy,t)w=logspace(-1,1,100)';[m p]=bode(num,den,w);figure(6)subplot(211);semilogx(w,20*log10(m)); grid onsubplot(212);semilogx(w,p)grid on[gm,pm,wpc,wgc]=margin(sys)figure(7)margin(sys)figure(8)nyquist(sys)figure(9)nichols(sys)方法3:F kx x b x m =++•••125=++•••x x xx x x 4.02.02.0--=•••% 主程序wffc1.mt=[0 30];x0=[0;0];[tt,yy]=ode45(@myfun1,t,x0);figure(1)plot(tt,yy(:,1),':b',tt,yy(:,2),'-r') hold onplot(tt,0.2-0.2*yy(:,2)-0.4*yy(:,1),'-.k ')legend('位移','速度','加速度') title('微分方程的解')figure(2)plot(yy(:,1),yy(:,2))title('平面相轨迹')%写出函数文件myfun1.mfunction xdot=myfun1(t,x)xdot=[x(2);0.2-0.2*x(2)-0.4*x(1)];。
Matlab求解微分方程(组)及偏微分方程(组)
第四讲【2 】Matlab求解微分方程(组)理论介绍:Matlab求解微分方程(组)敕令求解实例:Matlab求解微分方程(组)实例现实运用问题经由过程数学建模所归纳得到的方程,绝大多半都是微分方程,真正能得到代数方程的机遇很少.另一方面,可以或许求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就请求我们必须研讨微分方程(组)的解法:解析解法和数值解法.一.相干函数.敕令及简介1.在Matlab中,用大写字母D表示导数,Dy表示y关于自变量的一阶导数,D2y表示y关于自变量的二阶导数,依此类推.函数dsolve用来解决常微分方程(组)的求解问题,挪用格局为:X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve用来解符号常微分方程.方程组,假如没有初始前提,则求出通解,假如有初始前提,则求出特解.留意,体系缺省的自变量为t2.函数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 用初始前提y 求解.(3)假如要获得微分方程问题在其他指准时光点012,,,,ft t t t 上的解,则令tspan012[,,,]f t t t t =(请求是单调的).(4)因为没有一种算法可以有用的解决所有的ODE 问题,为此,Matlab 供给了多种求解器solver,对于不同的ODE 问题,采用不同的solver.表1 Matlab 中文本文件读写函数解释:ode23.ode45是极其常用的用来求解非刚性的标准情势的一阶微分方程(组)的初值问题的解的Matlab 常用程序,个中:ode23采用龙格-库塔2阶算法,用3阶公式作误差估量来调节步长,具有低等的精度.ode45则采用龙格-库塔4阶算法,用5阶公式作误差估量来调节步长,具有中等的精度.3.在matlab敕令窗口.程序或函数中创建局部函数时,可用内联函数inline,inli ne函数情势相当于编写M函数文件,但不需编写M-文件就可以描写出某种数学关系.挪用inline函数,只能由一个matlab表达式构成,并且只能返回一个变量,不许可[u,v]这种向量情势.因而,任何请求逻辑运算或乘法运算以求得最终成果的场合,都不能运用inline函数,inline函数的一般情势为:FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如:(求解F(x)=x^2*cos(a*x)-b ,a,b是标量;x是向量)在敕令窗口输入: Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’);g= Fofx([pi/3 pi/3.5],4,1)体系输出为:g=-1.5483 -1.7259留意:因为运用内联对象函数inline不须要别的树立m文件,所有运用比较便利,别的在运用ode45函数的时刻,界说函数往往须要编辑一个m文件来单独界说,如许不便于治理文件,这里可以运用inline来界说函数.二.实例介绍1.几个可以直接用Matlab求微分方程准确解的实例例1 求解微分方程2 '2x y xy xe-+=程序:syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)例 2求微分方程'0xxy y e+-=在初始前提(1)2y e=下的特解并画出解函数的图形.程序:syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x’);ezplot(y)例 3求解微分方程组530tdx x y e dt dy x y dt ⎧++=⎪⎪⎨⎪--=⎪⎩在初始前提00|1,|0t t x y ====下的特解并画出解函数的图形.程序: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 auto2.用ode23.ode45等求解非刚性标准情势的一阶微分方程(组)的初值问题的数值解(近似解)例4求解微分方程初值问题2222(0)1dy y x x dx y ⎧=-++⎪⎨⎪=⎩的数值解,求解规模为区间[0,0.5].程序:fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-') 例5求解微分方程22'2(1)0,(0)1,(0)0d y dy y y y y dt dt μ--+===的解,并画出解的图形.剖析:这是一个二阶非线性方程,我们可以经由过程变换,将二阶方程化为一阶方程组求解.令12,,7dyx y x dt μ===,则121221212,(0)17(1),(0)0dx x x dtdx x x x x dt ⎧==⎪⎪⎨⎪=--=⎪⎩编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 敕令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)演习与思虑:M-文件vdp.m 改写成inline 函数程序? 3.用Euler 折线法求解Euler 折线法求解的根本思惟是将微分方程初值问题00(,)()dyf x y dx y x y ⎧=⎪⎨⎪=⎩化成一个代数(差分)方程,重要步骤是用差商()()y x h y x h +-替代微商dydx ,于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例6用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解(步长h 取0.4),求解规模为区间[0,2].剖析:本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序:>> 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-1y=y+h*subs(f,{'x','y'},{x,y});%subs,调换函数 x=x+h; szj=[szj;x,y];end >>szj>> plot(szj(:,1),szj(:,2)) 解释:调换函数subs 例如:输入subs(a+b,a,4)意思就是把a 用4调换掉落,返回4+b,也可以调换多个变量,例如:subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分离用字符alpha 调换a 和2调换b,返回 cos(alpha)+sin(2)特别解释:本问题可进一步运用四阶Runge-Kutta 法求解,Euler 折线法现实上就是一阶Runge-Kutta 法,Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩响应的Matlab 程序为:>> 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-1l1=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))演习与思虑:(1)ode45求解问题并比较差异.(2)运用Matlab 求微分方程(4)(3)''20y y y -+=的解. (3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解.(4)运用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解.提示:尽可能多的斟酌解法三.微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准情势的一阶显式微分方程(组)问题,是以在运用ODE 解算器之前,我们须要做的第一步,也是最重要的一步就是借助状况变量将微分方程(组)化成Matlab 可接收的标准情势.当然,假如ODEs 由一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍若何将它变换成一个一阶显式微分方程组.Step1将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高分列.情势为:()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x xy y y y ----⎧=⎨=⎩ Step 2 为每一阶微分式选择状况变量,最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y --++++========留意:ODEs 中所有是因变量的最高阶次之和就是须要的状况变量的个数,最高阶的微分式不须要给它状况变量.Step 3 依据选用的状况变量,写出所有状况变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m n m n x x x x x x x f t x x x x x x x g t x x x x +++++======演习与思虑:(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩个中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩提示:运用符号盘算函数solve 求'''',x y ,然后运用求解微分方程的办法四.偏微分方程解法Matlab 供给了两种办法解决PDE 问题,一是运用pdepe 函数,它可以求解一般的PDEs,具有较大的通用性,但只支撑敕令情势挪用;二是运用PDE 对象箱,可以求解特别PDE 问题,PDEtoll 有较大的局限性,比如只能求解二阶PDE 问题,并且不能解决片微分方程组,但是它供给了GUI 界面,从庞杂的编程中摆脱出来,同时还可以经由过程File —>Save As 直接生成M 代码.1.一般偏微分方程(组)的求解(1)Matlab 供给的pdepe 函数,可以直接求解一般偏微分方程(组),它的挪用格局为:sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描写函数,它必须换成标准情势:(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x -∂∂∂∂∂=+∂∂∂∂∂如许,PDE 就可以编写进口函数:[c,f,s]=pdefun(x,t,u,du),m,x,t 对应于式中相干参数,du 是u 的一阶导数,由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界前提描写函数,它必须化为情势:(,,)(,,).*(,,,)0up x t u q x t u f x t u x ∂==∂于是边值前提可以编写函数描写为:[pa,qa,pb,qb]=pdebc(x,t,u,du),个中a 表示下边界,b 表示上边界.@pdeic 是PDE 的初值前提,必须化为情势:00(,)u x t u =,故可以运用函数描写为:u0=pdeic(x)sol 是一个三维数组,sol(:,:,i)表示i u 的解,换句话说,k u 对应x(i)和t(j)时的解为s ol(i,j,k),经由过程sol,我们可以运用pdeval 函数直接盘算某个点的函数值.(2)实例解释求解偏微分2111222221220.024()0.17()u u F u u t x u u F u u t x ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩个中, 5.7311.46()x x F x e e -=-且知足初始前提12(,0)1,(,0)0u x u x ==及边界前提1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0u u t u t t x ∂===∂解:(1)对比给出的偏微分方程和pdepe 函数求解的标准情势,原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ 可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦%目的PDE 函数function [c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp))end(2)边界前提改写为:下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ %边界前提函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];end(3)初值前提改写为:1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值前提函数function u0=pdeic(x)u0=[1;0];end(4)编写主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);subplot(2,1,1)surf(x,t,sol(:,:,1))subplot(2,1,2)surf(x,t,sol(:,:,2))演习与思虑: This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x x π∂∂∂=∂∂∂This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t u u t e t x π-∂=+=∂2.PDEtool 求解偏微分方程(1)PDEtool (GUI )求解偏微分方程的一般步骤在Matlab 敕令窗口输入pdetool,回车,PDE 对象箱的图形用户界面(GUI)体系就启动了.从界说一个偏微分方程问题到完成解偏微分方程的定解,全部进程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω,经由过程公式把Matlab 体系供给的实体模子:矩形.圆.椭圆和多边形,组合起来,生成须要的平面区域.Step 2 “Boundary 模式”界说边界,声明不同边界段的边界前提.Step3 “PDE 模式”界说偏微分方程,肯定方程类型和方程系数c,a,f,d,依据具体情况,还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω,可以掌握主动生成网格的参数,对生成的网格进行多次细化,使网格朋分更细更合理.Step 5 “Solve 模式”解偏微分方程,对于椭圆型方程可以激活并掌握非线性自顺应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,设置初始边界前提后可以求出给准时刻t 的解;对于特点值问题,可以求出给定区间上的特点值.求解完成后,可以返回到Step 4,对网格进一步细化,进行再次求解.Step 6 “View 模式”盘算成果的可视化,可以经由过程设置体系供给的对话框,显示所求的解的表面图.网格图.等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行为画演示.(2)实例解释用法求解一个正方形区域上的特点值问题:12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩正方形区域为:11,1 1.x x -≤≤-≤≤(1)运用PDE 对象箱打开GUI 求解方程(2)进入Draw 模式,绘制一个矩形,然后双击矩形,在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并封闭对话框(3)进入Boundary 模式,边界前提采用Dirichlet 前提的默认值(4)进入PDE 模式,单击对象栏PDE 按钮,在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后封闭对话框(5)单击对象栏的 按钮,对正方形区域进行初始网格剖分,然后再对网格进一步细化剖分一次(6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特点值区域为[ -20,20](7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color.Height(3-D plot)和show mesh项,然后单击Done确认(8)单击对象栏的“=”按钮,开端求解。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
包含有快速变化的部分。
如何调用?
y=dsolve('e1,e2,...','c1,c2,...','v')
其中'e1,e2,...'为微分方程或微分方程组; 'c1,c2,...',是初始条件或边界条件; 'v'是独立变量,默认的独立变量是't'; y 返回解析解。如果没有初始条件,则求出通解,如
ode23t、ode23tb 函数;
odefun 是函数句柄;
tspan 微分定义区间;
y0
为初值行矩阵;
T
值是t序列(为列向量);
Y
值是微分方程的解Y在各点t的值(为列向量);
TE
表示事件发生时间,可缺省;
YE
表示事件解决时间,可缺省;
IE
表示事件消失时间,可缺省;
options 是求解参数设置,可以用odeset在计算前设定误
Topic: 如何使用MATLAB求 解常微分方程(组)
TMU_BME_2013
a.What ?
微分方程指描述未知函数的导数与自变 量之间的关系的方程。未知函数是一元函 数的微分方程称作常微分方程。未知函数 是多元函数的微分方程称作偏微分方程。
MATLAB(matrix&laboratory)意为矩 阵工厂(矩阵实验室).MATLAB是美国 MathWorks公司出品的商业数学软件,提 供高级技术计算语言和交互式环境,主要 包括MATLAB和Simulink两大部分。
果有初始条件,则求出特解。 用字符串表示常微分方程,自变量缺省时为t,导数用
D表示微分。y的2阶导数用D2y表示,依此类推。
如olver('odefun',tspan,y0,options)
其中solver为ode23、ode45、ode113、ode15s、ode23s、
Where ?
工程控制 航空航天
ODE
医学生理 金融分析
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
When ?
当对问题进行建模后,有常微分方程 需要求解时。
在生物建模中,经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
How ?
数 值 解
数值解?刚性方程?
y(1)=x(2);
y1
y2
y(2)= -t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
1000
求解程序关键步骤
[t,y]=ode45('odefun',[3.9 4.0],[2 8])
y
500
0
3.9
3.92
3.94
3.96
3.98
4
4.02
t
Examples
E.g.3 求解人体中不同形式的碘浓度的三房室模型。已知碘在三房室之间的 转换速率为:k21=0.84/d,k01=1.68/d,k32=0.01/d,k13=0.08/d, k03=0.02/d,f10=150g/d,f30=0g/d。初始条件为:x1(0)=81.2g, x2(0)=6821g,x3(0)=682g,v1=4L, v2=2L, v3=5L。仿真时间为30d。
差,输出参数,事件等,可缺省。
使用ODE?时如何编 写微分方程 ?
方式一:带额外参数,使用时需对参数进行赋值
function odefun(t,x,flag,R,L,C) %用flag说明R、L、C为变
量
xdot=zeros(2,1);
%表明其为列向量
xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);
输入命令: [x,y,z]=dsolve(‘Dx=2*x-3*y+3*z’,’Dy=4*x-5*y+3*z’,
‘Dz=4*x4*y+2*z’,’t’) x=simple(x); y=simple(y); z=simple(z); %化简结果
运行结果为: x=(c1-c2+c3+c2e -3t-c3e-3t)e2t y= - c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1 e-4t+c2e-4t+c1-c2+c3)e2t
E.g.3结果
C1(t)/(ug/L)
20.3 20.295
20.29 20.285
0
3410.6
体 内 无 机 碘 浓 度 C1(t)
5
10
15
20
25
30
t/d
甲 状 腺 中 有 机 碘 浓 度 C2(t)
C2(t)/(ug/L)
Examples
E.g.2 求解y''== - t*y + e^t*y' +3sin(2t)。 已知3.9<t<4.0,y'(0)=2,y''(0)=8
函数文件odefun.m的建立
程序执行结果
function y=odefun(t,x) y=zeros(2,1); %列向量
y' '=-t*y + et*y' +3sin2t 1500
数值解是采用如有限元方 法、 数值逼近方法、插值方 法等计算方法得到的解。只 能利用数值计算的结果, 而 不能随意给出自变量并求出 计算值。
当无法藉由微积分技巧求 得解析解时,这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下,常 微分方程往往只能求解出其
数值解。
在数学中,刚性方程是指一 个微分方程,其数值分析的解 只有在时间间隔很小时才会稳 定,只要时间间隔略大,其解 就会不稳定。
通过列写方程组,建立odefun.m文 件即方程组文件 function dC=odefun(t,C)
dC=[0.1*C(3)+37.5-2.52*C(1); 1.68*C(1)-0.01*C(2); 0.004*C(2)-0.1*C(3)];
end 程序关键部分列写如下: [t,C]=ode15s('fun1',[0,30],[20.3,341 0.5,136.4])
xdot(2)=1/C*x(1);
e方n式d 二:无额外参数
function dC=odefun(t,C)
dC=[ 0.1*C(3)+37.5-2.52*C(1); %直接书写列矩阵
1.68*C(1)-0.01*C(2);
0.004*C(2)-0.1*C(3)];
end
Examples
E.g.1 求下列微分方程组的通解 x'=2x-3y+3z y'=4x-5y+3z z'=4x-4y+2z