MATLAB第四章符号运算2-微分方程求解
04第四章Matlab微分方程-Matlab教程
实际应用时，与欧拉公式结合使用：
yi(0) yi hf ( xi , yi ) 1 ( k 1) h k yi 1 yi [ f ( xi , yi ) f ( xi 1, yi(1) )] k 0,1, 2, 2
k k （k 对于已给的精确度，当满足 yi(11) yi(1) 时， yi 1 yi 11） 取 ，
则微分方程变为一阶微分方程组：
1、建立m-文件vdp1000.m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
-0.5 -1 -1.5 -2 -2.5 0
y1 ' y2 2 y2 ' 1000(1 y1 ) y2 y1 y (0) 2, y (0) 0 1 2
y1 ' y2 y3 y2 ' y1 y3 y3 ' 0.51 y1 y2 y1 (0) 0, y2 (0) 1, y3 (0) 1
1、建立m-文件rigid.m如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
2 1.5 1 0.5 0
500
1000
1500
2000
2500
3000
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-') 3、结果如图
MATLAB实验四_求微分方程的解
参数说明
[T,Y] = solver(odefun,tspan,y0)
odefun 为显式常微分方程，可以用命令 inline 定义，或 在函数文件中定义，然后通过函数句柄调用。
dy 2 2 y 2 x 2x 求初值问题 的数值解，求解范 例： dx 围为 [0,0.5] y( 0 ) 1
dsolve的输出个数只能为一个 或 与方程个数相等。
只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。
Matlab函数数值求解
[T,Y] = solver(odefun,tspan,y0)
其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (列向量) 中返回的是分割点 的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解， 其列数等于因变量的个数。
数学实验
实验四
求微分方程的解
问题背景和实验目的
自牛顿发明微积分以来，微分方程在描述事物运 动规律上已发挥了重要的作用。实际应用问题通过 数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。 然而能够求得解析解的微分方程十分有限，绝大多 数微分方程需要利用数值方法来近似求解。 本实验主要研究如何用 Matlab 来计算微分方程 （组）的数值解，并重点介绍一个求解微分方程的 基本数值解法－－Euler折线法。
Runge-Kutta 方法
Euler 法与 R-K法误差比较
Matlab 解初值问题
用 Maltab自带函数 解初值问题 求解析解：dsolve 求数值解：
ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb
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教程4
4.4
ቤተ መጻሕፍቲ ባይዱ拟合与插值
说明：interp1仅是 插值指令的一种， 还有interp2 、 interp3等。
ys=interp1(x,y,xs,’method’); 其中:(1) x,y是测量数据对； (2) xs是需要内插的点所构成的向量。 (3) method是指所使用的内插方法。 插值算法： ‘nearest’，‘linear’，‘spline’，‘cubic’
f ( x) g ( x) f ( x) / g ( x)
杨惠--《matlab语言及应用》--第四章
4.3 多项式计算
3.
功能：按数组运算规则计算多项式的值。 多项式值的计算： 其中x可以是标量和数组。
y=polyval(p,x); y=polyvalm(p,x);
功能：按矩阵运算规 则计算多项式的值。 其中x必须为方阵。
conv(p1,p2)
3、多项式除法运算
[Q,r]=deconv(p1,p2)
p1=conv(p2,Q)+r
例题开讲
f ( x) 3x 5 5 x 4 2 x 3 7 x 2 5 x 6 g ( x) 3x 2 5 x 3
已知
求： (1) f ( x) g ( x) (2) f ( x) g ( x)
4.4 拟合与插值
1）最邻近插值方法（nearest）
插值点的值与其最邻近的点的函数值相等。
*
*
*
*
*
4.4 拟合与插值
2）线性插值方法（ linear ）
插值点的值在前，后两个数据点所构 成的直线上。
*
*
*
*
*
4.4
拟合与插值
MATLAB8.5教程第4章 符号计算
• MATLAB为符号计算提供了一种引入符号对象的数学运 算工具箱，包含函数的复合、简化、极限、导数、积分， 泰勒展开式、级数求和，以及求解代数方程和微分方程 等函数命令。其计算指令的调用比较简单，基本上与数 学函数表示法相同。
本章重点
• • • • 符号对象的创建 符号极限、导数、积分 方程求解 级数求和
4.1.2 符号表达式运算
• 8．复合函数的运算 • 格式：compose(f,g) %返回f=f(x)和g=g(y)的复合函数 f(g(y))。x是为findsym定义的f的符号变量，y是为 findsym定义的g的符号变量 • compose(f,g,t) %返回f=f(x)和g=g(y)的复合函数f(g(t))， 返回的函数以t为自变量。x是为findsym定义的f的符号变 量，y是为findsym定义的g的符号变量。例如， • >>syms x y t • >>f=1/(1+x^2) • >>g=sin(y) • >>compose(f,g) • ans = • 1/(sin(y)^2 + 1) • >> compose(f,g,t) • ans = • 1/(sin(t)^2 + 1)
4.5 符号级数
• 2．级数和
• 格式：S=symsum(f) %对符号表达式f中的符号变量k （由命令findsym(f)确定的）从0到k-1求级数和
• S=symsum(f,x) 0到k-1求级数和
%对符号表达式f中指定的符号变量x从
• S=symsum(f,a,b) %对符号表达式f中的符号变量k（由 命令findsym(f)确定的）从a到b求级数和 • S=symsum(f,x,a,b) %对符号表达式f中指定的符号变量x 从a 到b求级数和
matlab利用ode45求解二元二阶微分方程
题目：探究matlab利用ode45求解二元二阶微分方程的方法与应用在数学和工程领域，微分方程是一类重要的数学工具，它可以描述自然界中众多的现象和规律。
而求解微分方程的问题一直是科学家和工程师们所关注的重要问题之一。
在计算机辅助数学建模领域，matlab作为一种强大的数值计算工具，可以通过内置的函数ode45来求解常微分方程初值问题。
本文将探讨matlab利用ode45求解二元二阶微分方程的方法与应用。
一、二元二阶微分方程的基本概念二元二阶微分方程是指含有两个自变量、二阶导数和一阶导数的微分方程。
一般形式如下：\[ F(x, y, \frac{dy}{dx}, \frac{d^2y}{dx^2}) = 0 \]其中x为自变量，y为因变量，\(\frac{dy}{dx}\)为y关于x的一阶导数，\(\frac{d^2y}{dx^2}\)为y关于x的二阶导数。
二、matlab中ode45函数的基本原理在matlab中，ode45是求解常微分方程初值问题的函数，它使用了一种自适应步长的Runge-Kutta方法来求解微分方程。
ode45可以求解一阶或高阶的常微分方程组，是matlab中最常用的求解微分方程的函数之一。
对于二元二阶微分方程，可以通过一些简单的变换和处理，转化为一组一阶微分方程的形式，然后利用ode45进行求解。
三、matlab利用ode45求解二元二阶微分方程的具体步骤1. 将二元二阶微分方程转化为一组一阶微分方程。
对于形如\(\frac{d^2y}{dx^2} = f(x, y, \frac{dy}{dx})\)的二阶微分方程，可以引入新的变量z = \(\frac{dy}{dx}\)，转化为一组一阶微分方程：\[\frac{dy}{dx} = z\]\[\frac{dz}{dx} = f(x, y, z)\]2. 编写matlab脚本文件。
在matlab中，编写脚本文件来定义微分方程的函数形式，并调用ode45函数来求解微分方程。
4MATLAB符号计算
第四节MATLAB符号计算在自然科学的各个领域不但需要解决数值分析和计算问题，同时也要解决符号运算的问题，MA TLAB中的符号计算功能是以Maple V为基础开发的。
MATLAB的符号数学工具箱的主要功能包括：符号表达式的创建、符号矩阵的运算，符号表达式的化简和替换、符号微积分、符号代数方程等。
一、符号表达式的创建MATLAB的符号数学工具箱提供了两个基本函数，用来创建符号变量和表达式，分别是sym 和syms。
●函数sym的调用形式为：x=sym(‘x’)创建一个符号变量x，它可以是字符、字符串、表达式或字符表达式。
●函数syms用于方便地一次创建多个符号变量，其调用形式为：syms a b c…例1 使用sym 和syms函数创建符号变量。
a=sym('a') %定义符号变量aa =ab=sym('1+sqrt(5)/2') %定义符号变量bb =1+sqrt(5)/2syms a b c d %定义4个符号变量使用函数可以创建符号矩阵，可以直接输入或从数值矩阵转换。
例2 创建一个循环矩阵。
syms a b c dn=[a b c d ;b c d a ; c d a b ; d a b c]输出结果为：n =[ a, b, c, d][ b, c, d, a][ c, d, a, b][ d, a, b, c]例3 将3阶的Hilbert 矩阵转化为符号矩阵。
h=hilb(3) %创建Hilbert矩阵h =1.0000 0.5000 0.33330.5000 0.3333 0.25000.3333 0.2500 0.2000h1=sym(h) %用函数sym转化为符号矩阵h1 =[ 1, 1/2, 1/3][ 1/2, 1/3, 1/4][ 1/3, 1/4, 1/5]注意：符号矩阵与普通数值矩阵的区别是：在命令窗口的显示中，数值矩阵只显示元素的数值，而符号矩阵的每行元素均放在一对方括号内；在工作空间窗口显示的变量图标也不同，数值图标为，符号矩阵的图标为。
Matlab教学第四章 MATLAB符号运算(Symbolic)
>> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') >> syms x; diff(y)+2*x*y - x*exp(-x^2)
f2=2*(u+2)
ans=14 ans=2*((a+2)+2) f3=2*x+2*y ans=6
符号矩阵
使用 sym 函数直接生成 >> A=sym('[1+x, sin(x); 5, exp(x)]') 将数值矩阵转化成符号矩阵 >> B=[2/3, sqrt(2); 5.2, log(3)]; >> C=sym(B) 符号矩阵中元素的引用和修改 >> A=sym('[1+x, sin(x); 5, exp(x)]'); >> A(1,2) % 引用 >> A(2,2)=sym('cos(x)') % 重新赋值
符号对象的基本运算
基本函数
三角函数与反三角函数、指数函数、对数函数等
sin、cos、tan、cot、sec、csc、… asin、acos、atan、acot、asec、 acsc、…
exp、log、log2、log10、sqrt abs、conj、real、imag
rank、det、inv、eig、lu、qr、svd
How 中记录的为简化过程中使用的方法。
f
2*cos(x)^2sin(x)^2
(x+1)*x*(x-1)
R
HOW
3*cos(x)^2-1 simplify
x^3-x combine(tri g)
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()xx 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符号计算
Ezpolar函数实现极坐标中二维曲 线图象的绘制,它的调用格式:
ezpolar(f,[a,b])
说明: 绘图表达式为rho=f(theta)极坐标曲线,
theta的取值范围为,缺省时为[0,2*pi]
例:p83
4.5.2三维绘图函数
ezplot3是实现三维绘图的函数,它的 调用格式:
ezplot3(x,y,z,[tmin,tmax]) ezplot3(x,y,z,[tmin,tmax ],'animate')
3.符号表达式的提取分子和分母运算
在matlab中可利用numden函数来提取符号表达式 中的分子或分母。其一般调用格式为:
[n,d]=numden(s)
说明:
参数s是符号表达式是一个 有理分式或可以展开为有理
例: s=sym('2/5+3/7')
n= 29
分式,numden函数把有理分 式的分子返回给n,分母返回
dsolve在求常微分方程组时 的调用格式为:
dsolve('eq1,eq2,...','cond1,cond2...','v')
说明: 该函数求解常微分方程组eq1,…,eqn在初值条件或边界 条件为cond1,…,condn下的符号特解，若不给出初值 条件，则求方程组的通解，v为指定的返回值中的变量 表示。
例4-9解超越方程组
sin(x+y)-y*exp(x)=0 x^2-y=2
syms x y
[x,y]=solve('sin(x+y)-y*exp(x)=0','x^2-y=2')
x =-.66870120500236202933135901833637
matlab求解二阶微分方程代码
在MATLAB中求解二阶微分方程是一个常见且关键的问题，尤其是在工程和科学领域中。
MATLAB提供了强大的求解工具和函数，可以帮助我们高效地解决这类问题。
我们需要明确二阶微分方程的一般形式：\[ a\frac{{d^2y}}{{dt^2}} + b\frac{{dy}}{{dt}} + cy = f(t) \]其中，\( a \)、\( b \) 和 \( c \) 是常数，\( y \) 是未知函数，\( f(t) \) 是已知函数。
接下来，我们使用MATLAB中的符号计算工具箱来求解二阶微分方程。
我们需要定义未知函数 \( y \) 作为符号变量，然后使用符号计算工具箱中的 \texttt{dsolve} 函数来求解微分方程。
```matlabsyms y(t)eqn = a*diff(y, t, 2) + b*diff(y, t) + c*y == f(t);ySol(t) = dsolve(eqn);```在上面的代码中，我们首先定义了未知函数 \( y \) 作为符号变量。
我们将二阶微分方程表示为等式 \( eqn \)，并使用 \texttt{dsolve} 函数来求解微分方程，得到了 \( y(t) \) 的解析解。
除了使用符号计算工具箱来求解微分方程外，我们还可以使用数值方法来求解。
MATLAB中提供了丰富的数值求解函数，如\texttt{ode45}、\texttt{ode23} 等，可以有效地求解微分方程的数值解。
```matlab[t, y] = ode45(@(t,y) myODE(t,y), tspan, y0);```在上面的代码中，我们使用了 \texttt{ode45} 函数来求解微分方程的数值解。
其中，\texttt{@(t,y) myODE(t,y)} 是对微分方程进行数值求解的函数句柄，\texttt{tspan} 是时间范围，\texttt{y0} 是初始条件。
matlab微分方程组的解法
一、引言1.1 MATLAB在微分方程组求解中的应用MATLAB作为一种强大的数学工具，被广泛应用于微分方程组的求解与模拟分析。
1.2 本文的研究目的和意义本文旨在探讨MATLAB在求解微分方程组方面的应用方法，帮助读者更好地理解和运用MATLAB进行微分方程组的解法，从而提高数学建模和工程仿真的效率与精度。
二、微分方程组的基本概念2.1 微分方程组的定义微分方程组是由多个未知函数及其偏导数构成的方程组。
常见的微分方程组可以分为线性微分方程组与非线性微分方程组。
2.2 微分方程组的求解方法求解微分方程组的方法包括解析解法、数值解法和符号解法。
而MATLAB在微分方程数值解法中具有独特的优势。
三、MATLAB在微分方程组求解中的基本操作3.1 MATLAB中微分方程组的表示在MATLAB中，微分方程组可以使用符号表达式或者函数形式表示，便于进行数值求解和仿真分析。
3.2 MATLAB中微分方程组的数值求解利用MATLAB中的ode45、ode23等求解微分方程组的函数，可以快速地求得微分方程组的数值解，并且可以灵活地控制求解的精度和速度。
3.3 MATLAB中微分方程组的图像绘制MATLAB提供了丰富的绘图函数，能够直观地展现微分方程组的数值解，帮助用户更直观地理解微分方程组的解法结果。
四、 MATLAB在微分方程组求解中的应用实例4.1 简单的线性微分方程组求解通过一个简单的线性微分方程组的求解实例，展示MATLAB在微分方程组求解中的基本操作和方法。
4.2 复杂的非线性微分方程组求解通过一个包含非线性项的微分方程组求解实例，展示MATLAB在处理复杂微分方程组时的应用能力。
五、MATLAB在微分方程组求解中的进阶应用5.1 高阶微分方程组的数值求解MATLAB可以利用符号运算工具箱对高阶微分方程组进行符号求解，也可以通过数值求解的方式得到高阶微分方程组的数值解。
5.2 特定约束条件下的微分方程组求解MATLAB可以通过引入特定的约束条件，对微分方程组进行求解，满足实际应用中的各种约束条件。
用MATLAB解常微分方程
实验四求微分方程的解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程， 真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限 的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组） 的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组） 的数值解法（近似解）.对微分方程（组）的解析解法（精确解），Matlab 有专门的函数可以用，本实 验将作一定的介绍.本实验将主要研究微分方程（组）的数值解法（近似解），重点介绍Euler 折线 法.二、相关函数（命令）及简介1. dsolve （'equ1','equ2'；…）：Matlab 求微分方程的解析解.equ1、equ2、… 为方程（或条件）.写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数， 用用D2y 表示y 关于自变量的二阶导数，依此类推.2. simplify （s ）:对表达式s 使用maple 的化简规则进行化简.例如：syms xsim pl ify （si n （x ）A2 + cos （x ）^2） an s=13. [r,how]=simple （s ）:由于 Matlab 提供了多种化简规则，simpie 命令就是 对表达式s 用各种规则进行化简，然后用r 返回最简形式，how 返回形成这种 形式所用的规则.例如：syms x[r,how]=sim ple （cos （x ）A2-si n （x ）A2） r = cos （2*x ） how = comb ine4. [T,Y] = solver （ odefun,tspan,y ））求微分方程的数值解.说明：（1）其中的 solver 为命令 ode45、ode23、ode113 ode15s ode23s 、ode23t 、 ode23tb 之一.⑶ 在积分区间tspan=[t 0,t f ]上，从t o 到t f ，用初始条件y 求解.⑵odefu n 是显式常微分方程:詈 f(t ，y)y(t 0) y o⑷ 要获得问题在其他指定时间点t0,t1,t2,上的解，则令tspan=[tott, ,t f ](要求是单调的).(5)因为没有一种算法可以有效地解决所有的ODE问题，为此，Matlab提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver.⑹要特别的是：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的Matlab的常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度.ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，有中等的精度.5. ezplot(x,y,[tmin,tmax]):符号函数的作图命令.x,y为关于参数t 号函数，[tmin,tmax]为t的取值范围.6.iniine():建立一个内联函数.格式：inline('expr', 'var1', 'var2',…)，括号里的表达式要加引号.例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)的符注意x三、实验内容1.几个可以直接用Matlab 求微分方程精确解的例子：例1:求解微分方程dy2xy xe dx求解本问题的Matlab 程序为：syms x yy=dsolve('Dy+2*x*y=x*ex p(-xA2)','x') diff(y,x)+2*x*y-x*ex 卩(帜八2)sim plify(diff(y,x)+2*x*y-x*ex 卩(帜八2))说明：(1)行Iine1是用命令定义x,y 为符号变量.这里可以不写,但为确保正确性, 建议写上；⑵ 行Iine2是用命令求出的微分方程的解：1/2*ex p(-xA2)*xA2+ex p(-xA2)*C1(3)行Iine3使用所求得的解.这里是将解代入原微分方程，结果应该为0, 但这里给出：-x^3*ex p(-xA2)-2*x*ex p(-xA2)*C1+2*x*(1/2*ex p(-xA2)*xA2+ex p(-xA2)*C1)确是微分方程的解.e x 0在初始条件y(1) 2e 下的特解，并画出解函 数的图形.求解本问题的Matlab 程序为：syms x yy=dsolve('x*Dy+y-ex p(x)=0','y(1)=2*ex p(1)', 'x')ezpl ot(y)e e x微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式),即卩 y --------- 解函数的图形如图1:2 x，并加以验证.%li ne1 %li ne2 %li ne3 %li ne4⑷行 Iine4 用 simplify()函数对上式进行化简，结果为0,表明y y(x)的例2:求微分方程xy' y并画出解函数的图形.求解本问题的Matlab 程序为：syms x y t[x,y]=dsolve('Dx+5*x+y=ex p(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') sim ple(x); sim ple(y);ezpl ot(x,y,[0,1.3]);axis auto微分方程的特解(式子特别长)以及解函数的图形均略.2. 用ode23 ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问 题的数值解(近似解)•例4：求解微分方程初值问题dy22y(0) 1间[0, 0.5].fun=i niin eC-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fu n, [0,0.5],1); x'; y';plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.2900 0.34000.39000.4400>> y'-30dx例3:求微分方程组dtdy5xx 3ye t在初始条件x|t 0 1,yl t 0 0下的特解,2x 2x的数值解，求解范围为区0.1900 0.2400 0.49000.50001/x exp(x)+1/x exp(1)5040302010-10-20-6-4-2246x0.9247 0.8434 0.7754 0.7199 0.6764 0.62220.61050.6084 0.6154 0.61792.ans =1.0000 0.6440图形结果为图0.6例5:求解描述振荡器的经典的Ver der Pol 微分方程分析：令 d 2y dt 2y 0,y(0) 1, y'(0) 0, 7.X 1 y ,x 2dx 1 则 dx 1 dt ， dtdx 2 X 2,~dr2(1 X 1 )X 2 X 1.先编写函数文件verderpol.m ： fun ctio nxp rime = verder pol(t,x) global mu; xp rime = [x(2) ;mu*(1-x(1)^2)*x (2) -x(1)]; 再编写命令文件vdp 1.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.0.950.90.850.80.750.70.650.150.20.250.30.350.40.450.50.053. 用Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的.下面介绍用Euler 折线 法求微分方程的数值解（近似解）的方法.Euler 折线法求解的基本思想是将微分方程初值问题dx S y(x 0) y 。
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符号微积分微分⽅程符号极限、微积分和符号⽅程的求解1.语法：sym(‘表达式’)%创建符号表达式f1=sym('a*x^2+b*x+c')f1 =a*x^2+b*x+c2.使⽤syms命令创建符号变量和符号表达式语法：syms arg1 arg2 …,参数%把字符变量定义为符号变量的简洁形式syms a b c x %创建多个符号变量f2=a*x^2+b*x+c %创建符号表达式3.4.1符号极限假定符号表达式的极限存在，Symbolic Math Toolbox提供了直接求表达式极限的函数limit，函数limit的基本⽤法如表3.2所⽰。
【例3.14】分别求1/x在0处从两边趋近、从左边趋近和从右边趋近的三个极限值。
f=sym('1/x')limit(f,'x',0) %对x求趋近于0的极限ans =NaNlimit(f,'x',0,'left') %左趋近于0ans =-inflimit(f,'x',0,'right') %右趋近于0ans =inf程序分析：当左右极限不相等，表达式的极限不存在为NaN。
3.4.2符号微分函数diff是⽤来求符号表达式的微分。
语法：diff(f) %求f对⾃由变量的⼀阶微分diff(f,t) %求f对符号变量t的⼀阶微分diff(f,n) %求f对⾃由变量的n阶微分diff(f,t,n) %求f对符号变量t的n阶微分【例3.15】已知f(x)＝ax2+bx+c，求f(x)的微分。
f=sym('a*x^2+b*x+c')f =a*x^2+b*x+cdiff(f) %对默认⾃由变量x求⼀阶微分ans =2*a*x+bdiff(f,'a') %对符号变量a求⼀阶微分ans =x^2diff(f,'x',2) %对符号变量x求⼆阶微分ans =2*adiff(f,3) %对默认⾃由变量x求三阶微分ans = 0微分函数diff 也可以⽤于符号矩阵，其结果是对矩阵的每⼀个元素进⾏微分运算。
matlab算法-求解微分方程数值解和解析解
MATLAB是一种用于数学计算、工程和科学应用程序开发的高级技术计算语言和交互式环境。
它被广泛应用于各种领域，尤其在工程和科学领域中被用于解决复杂的数学问题。
微分方程是许多工程和科学问题的基本数学描述，求解微分方程的数值解和解析解是MATLAB算法的一个重要应用。
1. 求解微分方程数值解在MATLAB中，可以使用各种数值方法来求解微分方程的数值解。
其中，常见的方法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。
这些数值方法可以通过编写MATLAB脚本来实现，从而得到微分方程的近似数值解。
以常微分方程为例，可以使用ode45函数来求解微分方程的数值解。
该函数是MATLAB中用于求解常微分方程初值问题的快速、鲁棒的数值方法，可以有效地得到微分方程的数值解。
2. 求解微分方程解析解除了求解微分方程的数值解外，MATLAB还可以用于求解微分方程的解析解。
对于一些特定类型的微分方程，可以使用符号计算工具箱中的函数来求解微分方程的解析解。
通过符号计算工具箱，可以对微分方程进行符号化处理，从而得到微分方程的解析解。
这对于研究微分方程的性质和特点非常有帮助，也有助于理论分析和验证数值解的准确性。
3. MATLAB算法应用举例在实际工程和科学应用中，MATLAB算法求解微分方程问题非常常见。
在控制系统设计中，经常需要对系统的动态特性进行分析和设计，这通常涉及到微分方程的建模和求解。
通过MATLAB算法，可以对系统的微分方程进行数值求解，从而得到系统的响应曲线和动态特性。
另外，在物理学、生物学、经济学等领域的建模和仿真中，也经常需要用到MATLAB算法来求解微分方程问题。
4. MATLAB算法优势相比于其他数学软件和编程语言，MATLAB在求解微分方程问题上具有明显的优势。
MATLAB提供了丰富的数值方法和工具，能够方便地对各种微分方程进行数值求解。
MATLAB具有直观的交互式界面和强大的绘图功能，能够直观地展示微分方程的数值解和解析解，有利于分析和理解问题。
Matlab教学第四章MATLAB符号运算(Symbolic)
符号计算可以给出完全正确的封闭解，或任意精度的数
值解（封闭解不存在时）。
符号计算指令的调用比较简单，与数学教科书上的公式相近。
符号计算所需的运行时间相对较长。
Matlab 符号运算举例
求一元二次方程 ax2 + bx + c = 0 的根 >> solve('a*x^2+b*x+c=0') 求的根 f (x) = (cos x)2 的一次导数 >> x=sym('x'); >> diff(cos(x)^2) 计算 f (x) = x2 在区间 [a, b] 上的定积分
>> syms a b x; >> int(x^2,a,b)
符号对象与符号表达式
在进行符号运算时，必须先定义基本的符号对象，可以是符号常量、符号变 量、符号表达式等。符号对象是一种数据结构。
含有符号对象的表达式称为符号表达式，Matlab 在内部把符号表达式表示 成字符串，以与数字变量或运算相区别。
syms a b c
符号表达式的建立
符号表达式的建立：
建立符号表达式通常有以下2种方法： (1) 用 sym 函数直接建立符号表达式。 (2) 使用已经定义的符号变量组成符号表达式。
例： >>
y=sym('sin(x)+cos(x)')
>> x=sym('x'); >> y=sin(x)+cos(x) >> syms x; >> y=sin(x)+cos(x)
六类常见符号运算
因式分解、展开、合并、简化及通分等 计算极限 计算导数 计算积分 符号求和 代数方程和微分方程求解
matlab使用欧拉法求解二阶微分代码
欧拉法是一种简单的数值方法，用于求解微分方程。
以下是使用MATLAB求解二阶微分方程的欧拉法代码：
```matlab
定义参数和初始条件
a = 0; 区间左端点
b = 1; 区间右端点
h = 0.01; 步长
x = a:h:b; 生成离散点
y = zeros(size(x)); 初始化y数组
y(1) = 1; 初始条件
y(2) = 0;
使用欧拉法进行迭代
for i = 2:length(x)
y(i) = y(i-1) + h * f(x(i-1), y(i-1)); f 是二阶微分方程的右端函数
end
定义微分方程的右端函数
f = @(x, y) x^2 - y;
绘制结果图
plot(x, y);
xlabel('x');
ylabel('y');
```
注意，这个代码只适用于简单的二阶微分方程，例如`y'' = f(x, y)`。
对于更复杂的微分方程，你可能需要修改`f` 函数以适应你的需求。
同时，这个代码没有包含任何错误处理或边界检查，所以在实际使用时可能需要添加这些内容。
matlab符号微积分微分方程
matlab符号微积分微分方程符号极限、微积分和符号方程的求解 1.语法:sym(…表达式?) %创建符号表达式f1=sym('a*x^2+b*x+c')f1 =a*x^2+b*x+c2.使用syms命令创建符号变量和符号表达式语法:syms arg1 arg2 …,参数 %把字符变量定义为符号变量的简洁形式syms a b c x %创建多个符号变量f2=a*x^2+b*x+c %创建符号表达式3.4.1符号极限假定符号表达式的极限存在，Symbolic Math Toolbox提供了直接求表达式极限的函数limit，函数limit的基本用法如表3.2所示。
表3.2 limit函数的用法表达式函数格式说明limf(x) limt(f) 对x求趋近于0的极限 x,0limf(x) limt(f,x,a) 对x求趋近于a的极限，当左右极限不相等时x,a极限不存在。
limf(x) limt(f,x,a, 对x求左趋近于a的极限 ,x,aleft)limf(x) limt(f,x,a, 对x求右趋近于a的极限，x,aright)【例3.14】分别求1/x在0处从两边趋近、从左边趋近和从右边趋近的三个极限值。
f=sym('1/x')limit(f,'x',0) %对x求趋近于0的极限ans =NaNlimit(f,'x',0,'left') %左趋近于0ans =-inflimit(f,'x',0,'right') %右趋近于0ans =inf程序分析:当左右极限不相等，表达式的极限不存在为NaN。
3.4.2符号微分函数diff是用来求符号表达式的微分。
语法:diff(f) %求f对自由变量的一阶微分diff(f,t) %求f对符号变量t的一阶微分diff(f,n) %求f对自由变量的n阶微分diff(f,t,n) %求f对符号变量t的n阶微分2【例3.15】已知f(x),ax+bx+c，求f(x)的微分。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
4/1
线性方程组的解析解
符号积分变换
求解一元二次方程a * x 2 + b * x + c = 0的解。 ¿¿ f=sym(’a*x^2+b*x+c=0’); ¿¿ xf=solve(f)
彭朝霞 (北京航空航天大学)
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
11 / 1
微分方程解析解
微分方程的解析解dsolve输出顺序
求微分方程组 du =v dt dv =w dt dw = −u dt
并且初值满足： u (0) = 0, v (0) = 0, w (0) = 1. 通解的MATLAB程序为： ¿¿S1 = dsolve(’Du=v, Dv=w, Dw=-u’,’u(0)=0, v(0)=0, w(0)=1’) ¿¿S2 = dsolve(’Dw=-u,Du=v, Dv=w’, ’w(0)=1, u(0)=0, v(0)=0’) ¿¿S3 = dsolve(’Dw=-u,Dv=w, Du=v’, ’w(0)=1, v(0)=0, u(0)=0’) 注： 解微分方程组时， 如果所给的输出个数与方程个数相同， 则方程组的解按词典 顺序输出；
彭朝霞 (北京航空航天大学)
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
6/1
线性方程组的解析解
求下列线性代数方程组的解。
x1 cos(sita) − x2 sin(sita) = a x1 sin(sita) + x2 cos(sita) = b
¿¿ syms x1 x2 a b sita; ¿¿ L1=x1 * cos(sita)-x2 * sin(sita)-a; ¿¿ L2=x1 * sin(sita)+x2 * cos(sita)-b; ¿¿ [x1,x2]=solve(L1,L2,x1,x2) 运行结果： x1 = cos(sita)*a+b*sin(sita) x2 = -a*sin(sita)+cos(sita)*b
1
微分方程中用D 表示对自变量的导数/微分， 如： 求高阶微分： Dy → y ′ ; D 2y → y ′′ ; D 3y → y ′′′ 。任何D 后面所跟的字母为因变量。
2 3
如果省略初值条件， 则表示求通解； 如果省略自变量， 则默认自变量为t dsolve(’Dy=2*x’,’x’); % dy/dx = 2x dsolve(’Dy=2*x’); % dy/dt = 2x
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
8/1Biblioteka 微分方程解析解微分方程的解析解
函数dsolve用来解符号微分方程、 方程组。 dsolve 的使用 y=dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
如果没有初始条件， 则求出通解， 如果有初始条件， 则求出特解。
14 / 1
微分方程数值解
数值解法的基本思想
xi +1 − xi = h, i = 0, 1, 2, . . . , n − 1,用离散的方法求解微分方程： y ˙ = f (x , y ) y (x0 ) = y0 ♣ 欧拉方法,即用差商代替导数 ♣ 改进的欧拉方法 ♣ 使用泰勒公式 ♣ 泰勒公式法： 龙格-塔库法、 线性多步法等方法
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
3/1
线性方程组的解析解
代数方程的求解solve
解方程x 3 − 3 * x + 1 = 0. f 是符号表达式 ¿¿ syms x; f=x^3-3*x+1; s=solve(f,x) 注意： 不能写成s=solve(f=0,x) f 是字符串 ¿¿ s=solve(’x^3-3*x+1’,’x’) ¿¿ s=solve(’x^3-3*x+1=0’,’x’) 注意： 可以不含等号， 表示解方程f=0， 但是也可以含等号。
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
13 / 1
微分方程数值解
2、 微分方程的数值解
♣ 常微分方程数值解的定义: 在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解.而在实际上对 初值问题， 一般是要求得到解在若干个点上满足规定精确度的近似值， 或者得到一 个满足精确度要求的便于计算的表达式. 因此， 研究常微分方程的数值解法是十分必要的. 对常微分方程 y ˙ = f (x , y ), y (x0 ) = y0
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
17 / 1
微分方程数值解
微分方程求数值解
第二步： 解微分方程 tfinal=100; x0=[0 0 1e-10]; [T , X ]=ode23(’zmfun’,[0,tfinal],x0)； plot(t,x)
MATLAB基础及工程应用
2014 年 10 月 23 日
5/1
线性方程组的解析解
代数方程求解
solve 也可以用来解方程组 solve (’eqn1’,’eqn2’,…,’eqnN’,’var1,var2,…,varN’) 求解由’eqn1’,’eqn2’,…,’eqnN’ 组成的方程组关于变量’var1,var2,…,varN’ 的解例 如： ¿¿ [x,y,z]=solve(’x+2*y-z=27’,’x+z=3’, ... ’x^2+3*y^2=28’,’x’,’y’,’z’) 输出变量的顺序要书写正确！ solve 在得不到解析解时， 会给出数值解
MATLAB基础及工程应用
2014 年 10 月 23 日
2/1
线性方程组的解析解
线性方程组的解析解
♣ 代数方程（非线性）求解： solve
s=solve(f,v) ： 求方程关于指定自变量的解 s=solve(f) ： 求方程关于默认自变量的解
f 可以是用字符串， 或符号表达式 来表示的方程
若f 是字符串， 可以不含等号， 表示解方程f=0 若f 是符号表达式， 不能含等号
[X,Y]=dsolve(’Dx+5*x+y=exp(t)’,’Dy-x-3*y=0’, ... ’x(0)=1’, ’y(0)=0’, ’t’) X =exp(2*t)/6 + (2*exp(t))/11 - 4*C3*exp(15^(1/2)*t - t) (4*C4)/exp(t + 15^(1/2)*t) + 15^(1/2)*C3*exp(15^(1/2)*t - t) (15^(1/2)*C4)/exp(t + 15^(1/2)*t) Y =C3*exp(15^(1/2)*t - t) - exp(t)/11 - (7*exp(2*t))/6 + C4/exp(t + 15^(1/2)*t) 注： 解微分方程组时， 如果所给的输出个数与方程个数相同， 则方程组的解按词典 顺序输出； 如果只给一个输出， 则输出的是一个包含解的结构（structure）类型的 数据。
第四章MATLAB 符号运算
彭朝霞 Email: pengzhaoxia@ 北京航空航天大学 交通科学与工程学院 2014 年 10 月 23 日
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
1/1
本章主要内容
彭朝霞 (北京航空航天大学)
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
15 / 1
微分方程数值解
2、 微分方程的数值解
♣ 微分方程数值求解函数与基本应用格式为： [t,y]=ode23(‘fun’,T0,Tfinal,Y0) [t,y]=ode45(‘fun’,T0,Tfinal,Y0) 函数说明： (1) ode23采用龙格-库塔二阶算法 ode45采用龙格-库塔四阶算法 (2) fun为由m函数定义的微分方程组的m函数名。 (3) T0,Tfinal,为数值积分区间， 仿真范围 (4) Y0为微分方程的初始条件向量。
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
7/1
线性方程组的解析解
求下列非线性代数方程组以y、 z 作变量的解。
uy 2 + vz + w = 0 y +z +w =0 ¿¿ syms y z u v w; ¿¿ eq1=u*y^2+v*z+w; ¿¿ eq2=y+z+w; ¿¿ [y z]=solve(eq1,eq2, y, z) y = [-1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))-w] [-1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))-w] z = [1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))] [1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))]
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
16 / 1
微分方程数值解
微分方程求数值解
例： 解下列微分方程组： 8 x ˙ 1 = − x1 + x2 x3 3 x ˙ 2 = −10x2 + 10x3 x ˙ 3 = −x1 x2 + 28x2 − x3 初始条件： x1 (0) = x2 (0) = 0, x3 = ������ = 1−10 。 1、 求微分方程的数值解 第一步： 建立M函数zmfum.m function dx=zmfun(t,x) dx=zeros(3,1); dx(1)=-8*x(1)/3+x(2)*x(3); dx(2)= -10*x(2)+10*x(3); dx(3)= -x(1)*x(2)+28*x(2)-x(3);