第四讲 MATLAB解方程及微分方程
MATLAB差异方程与微分方程求解技巧
MATLAB差异方程与微分方程求解技巧差异方程和微分方程是数学中重要的概念和工具，它们在各个领域的建模和分析中发挥着重要作用。
而MATLAB作为一款强大的数学软件，提供了丰富的工具和函数来求解差异方程和微分方程。
本文将介绍MATLAB中差异方程和微分方程的求解技巧，并提供一些实际案例来加深理解。
一、差异方程的求解技巧差异方程是描述离散域系统的数学模型，通常用递归关系来表达。
MATLAB 提供了多种方法来求解差异方程，其中最常用的是通过递推关系进行迭代。
1. 递推法递推法是通过迭代计算差异方程中的每一项来求解整个方程。
首先，需要定义差异方程的初始条件和递推关系。
然后，可以使用循环结构来进行迭代计算，直到达到所需精度或迭代次数。
假设我们要求解以下差异方程：y[n] = a * y[n-1] + b * y[n-2]其中，a和b为常数，y[n]为求解的项，y[n-1]和y[n-2]为已知的前两项。
在MATLAB中，可以使用for循环或while循环来实现递推法求解差异方程。
以下是使用for循环的实例代码：``` MATLABn = 1:10; % 定义计算的范围y = zeros(size(n)); % 初始化y的空间y(1) = y0; % 设定初始条件y(2) = y1; % 设定初始条件for i = 3:length(n)y(i) = a * y(i-1) + b * y(i-2); % 递推计算end```2. 齐次差异方程和非齐次差异方程的求解在求解差异方程时，需要区分齐次差异方程和非齐次差异方程。
对于齐次差异方程，它的非零解为零解；对于非齐次差异方程，它的非零解可以通过叠加齐次解和特解来得到。
MATLAB中，可以使用dsolve函数来求解差异方程。
以下是求解一阶齐次差异方程的实例代码：``` MATLABsyms y(t); % 定义符号变量eqn = diff(y, t) == a * y; % 定义差异方程cond = y(0) == y0; % 定义初始条件ySol(t) = dsolve(eqn, cond); % 求解差异方程```二、微分方程的求解技巧微分方程是描述连续域系统的数学模型，通常用导数关系来表达。
matlab求解微分代数方程
文章主题：探索数学求解软件Matlab在微分代数方程求解中的应用1. 引言微分代数方程（DAE）是描述物理系统中的相互依赖性和复杂性的数学模型。
解决这类方程对于现代科学和工程领域至关重要。
Matlab作为一种强大的数学计算软件，在微分代数方程求解中具有独特的优势。
本文将从简单到复杂的方式，探讨Matlab在DAE求解中的应用，并共享个人见解。
2. DAE的基本概念微分代数方程是描述包含未知函数及其导数或导数与未知函数的组合的方程。
通常的形式为F(x, x', t) = 0，其中x为未知函数，x'为其导数，t为自变量。
在实际应用中，这些方程往往伴随着初始条件和边界条件。
3. Matlab在解常微分方程（ODE）中的应用Matlab拥有丰富的ODE求解函数，如ode45、ode23等，可用于求解各种常微分方程。
这些函数可以自动选择适当的数值积分方法，并提供了方便的接口和参数设置，极大地简化了求解过程。
4. DAE求解方法的挑战与ODE相比，DAE的求解更具挑战性。
由于包含了代数变量和微分变量，常规的数值积分方法难以直接应用。
而且，方程的初始条件和边界条件也增加了求解的复杂性。
5. Matlab在DAE求解中的工具Matlab提供了一系列专门用于DAE求解的函数和工具包，如dare和ddesd等。
这些工具在模型建立、数值解法选择、收敛性分析等方面都具有独特的优势。
6. 案例分析：用Matlab求解电路模型的DAE以电路模型的DAE为例，通过Matlab可以快速建立系统的数学模型，并进行数值求解。
通过对参数的调节和模型的分析，可以更好地理解电路的动态特性，帮助优化设计和故障诊断。
7. 总结与展望通过本文的探讨，我们更深入地了解了Matlab在微分代数方程求解中的重要性和应用。
在未来，随着科学技术的发展，Matlab在此领域的功能和性能将得到进一步的提升，为工程科学领域提供更强大的支持。
个人观点：Matlab作为一种综合性的科学计算软件，对微分代数方程的求解起着至关重要的作用。
matlab求微分方程数值解
matlab求微分方程数值解利用matlab求微分方程数值解是一种常用的数学计算方法。
在实际工程和科学研究中，许多问题都可以用微分方程来描述，但是解析解往往难以求得，因此需要用数值方法求解微分方程。
求解微分方程的数值方法有很多种，其中比较常用的是欧拉法和龙格-库塔法。
欧拉法是一种基本的数值方法，它采用离散化的方法将微分方程转化为差分方程，然后通过迭代来求出数值解。
欧拉法的具体步骤是：首先将自变量和因变量离散化，然后利用微分方程的定义式将微分方程转化为差分方程，最后通过迭代求出数值解。
欧拉法的优点是简单易懂，但是精度较低，容易产生误差。
龙格-库塔法是一种高阶数值方法，它将微分方程转化为一系列的差分方程，并采用递推的方法求解数值解。
龙格-库塔法的优点是精度高，收敛速度快，适用于求解复杂的微分方程。
但是龙格-库塔法的计算量较大，需要进行多次计算，计算时间较长。
在使用matlab求解微分方程时，可以直接调用matlab中的ode 函数来求解微分方程。
ode函数是matlab中内置的求解微分方程的函数，它支持多种数值方法，包括欧拉法和龙格-库塔法等。
使用ode函数可以简化求解微分方程的过程，提高计算效率。
在使用ode函数求解微分方程时，需要先定义微分方程的函数表达式，然后将函数表达式作为参数传入ode函数中。
ode函数会自动选择合适的数值方法来求解微分方程，并返回数值解。
通过调整ode函数的参数，可以进一步提高求解微分方程的精度和计算效率。
除了ode函数外，matlab中还有很多其他的数值计算函数，如dsolve函数、pdepe函数等，它们可以用来求解不同类型的微分方程。
在实际应用中，需要根据具体问题选择合适的数值方法和函数来求解微分方程。
利用matlab求解微分方程数值解是一种常用的数学计算方法，可以通过调用matlab中的内置函数来实现。
在选择数值方法和函数时需要考虑精度和计算效率等因素，以便更好地解决实际问题。
matlab求解微分方程
matlab求解微分方程
Matlab是一种广受欢迎的工程数学软件，它可以解决复
杂的数学问题，其中包括求解常微分方程（ODEs）。
ODEs
在几何、物理和工程上都有很多应用，包括解决振动问题、热传递问题、控制系统设计问题等。
Matlab有一个内置的ODE解算器，可以根据给定的微分
方程求解。
Matlab中的ODE解算器可以分为两类：一类是ODE
45，它使用四阶Runge-Kutta法；另一类是ODE15s，它
使用单步积分器，如Adams-Bashforth和Adams-Moulton方法。
当使用Matlab求解ODE时，首先必须编写一个函数来定
义微分方程右端的函数，然后调用ODE解算器来求解，可以
指定边界条件、初始值、解算器类型等参数。
Matlab的可视化工具可以用来可视化求解的解，可以方
便地查看求解的结果，从而更好地理解问题。
总之，Matlab是一款功能强大的工程数学软件，它可以
很方便地求解ODE，可以用于各种科学和工程应用，是一款
非常强大的软件。
matlab解微分方程的通解
matlab解微分方程的通解
一、MATLAB解微分方程的通解
MATLAB有两种不同的函数可以解微分方程，一种是dsolve，另一种是ode45。
1.dsolve函数
dsolve函数是最常用的求解微分方程的函数，它可以求解一阶、二阶及更高阶的常微分方程，它能够得到方程的通解，但是它只能解指定类型的非线性方程，例如：常微分方程的通解，初值问题的解。
使用dsolve函数时应该按照以下步骤：
（1）输入微分方程，将其表示为一个字符串；
（2）调用dsolve函数，并传入字符串；
（3）调用结果，观察输出；
（4）如果输出为“未定义”，则需要检查输入的字符串是否正确；
（5）如果输出正确，则可以根据输出解析得到通解。
例1：解以下微分方程：
y″+2y′+5y=0
解：
首先将微分方程表示为字符串：
syms y
eqn=diff(y,2)+2*diff(y,1)+5*y==0
然后调用dsolve函数：
ySol=dsolve(eqn)
解析输出：
ySol=C1*exp(-3*t)+C2*exp(-2*t)
得到方程的通解为：
y=C1*exp(-3*t)+C2*exp(-2*t)
2.ode45函数
ode45函数是MATLAB中用于求解微分方程的另一种函数，它可以求解一阶、二阶及更高阶的常微分方程，以及积分方程、常系数线性微分方程等。
使用ode45函数时应该按照以下步骤：
（1）创建微分方程的函数；
（2）定义起始点和终止点；
（3）调用ode45函数，并传入函数及起始点和终止点；
（4）观察输出；
（5）根据结果获取通解。
matlab求微分方程组的解析解
MATLAB求微分方程组的解析解引言在科学与工程领域，微分方程组是一种常见的数学模型，用于描述各种物理现象和工程问题。
解析解是指能够用公式表达出来的精确解。
本文将介绍如何使用M ATL A B求解微分方程组的解析解。
1.方程组的建立首先，我们需要确定待求解的微分方程组。
假设我们有一个由n个微分方程组成的方程组，可以写为如下形式：d y1/dt=f1(t,y1,y2,...,yn)d y2/dt=f2(t,y1,y2,...,yn)......d y n/dt=f n(t,y1,y2,...,yn)其中`t`是自变量，`y1,y2,...,y n`是因变量，`f1,f2,...,fn`是给定的函数关系。
我们的目标是求解`y1(t),y2(t),...,yn(t)`的解析解。
2.使用MAT LAB求解M A TL AB提供了强大的求解微分方程组的工具，我们可以使用其中的函数来求解上述方程组的解析解。
首先，我们需要在MA T LA B中定义方程组的函数形式。
可以通过定义一个函数或者使用匿名函数来实现。
例如，我们可以定义一个名为`m yE qu at io ns`的函数，其输入参数为`t`和一个向量`y`，输出为一个向量`d y`，代表方程组的左侧和右侧的变量分别。
函数示例如下：f u nc ti on dy=m yE qua t io ns(t,y)%定义方程组d y=z er os(n,1);d y(1)=f1(t,y(1),y(2),...,y(n));d y(2)=f2(t,y(1),y(2),...,y(n));......d y(n)=fn(t,y(1),y(2),...,y(n));e n d然后，我们可以使用M AT LA B的`d so lv e`函数来求解微分方程组的解析解。
示例如下：s y ms ty1(t)y2(t)...yn(t)a s su me(t,'re al')a s su me(y1(t),'rea l')a s su me(y2(t),'rea l')......a s su me(y n(t),'rea l')e q n1=d if f(y1(t),t)==f1(t,y1(t),y2(t),...,y n(t));e q n2=d if f(y2(t),t)==f2(t,y1(t),y2(t),...,y n(t));......e q nn=d if f(yn(t),t)==fn(t,y1(t),y2(t),...,y n(t));e q ns=[eq n1,e qn2,...,eq nn];S=ds ol ve(e qn s);`S`即为方程组的解析解集合。
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求解微分方程特解。
二、微分方程特解的概念微分方程的一般形式可表示为：dy/dx = f(x, y)其中y是未知函数，x是自变量，f是已知函数。
微分方程的特解是指满足特定初值条件的解，通常表示为y(x0) = y0，其中x0和y0是已知的初值。
三、MATLAB求解微分方程特解的基本步骤1. 定义微分方程在MATLAB中，首先需要定义微分方程的函数形式。
假设我们要求解的微分方程为dy/dx = x + y，则在MATLAB中可以定义函数形式为：function dydx = myfun(x, y)dydx = x + y;2. 定义初值条件接下来需要定义初值条件，即给定的初始条件。
假设初值条件为y(0)= 1，则在MATLAB中可以定义为：x0 = 0;y0 = 1;3. 求解微分方程通过调用MATLAB中的内置函数ode45，可以求解微分方程的特解。
具体的求解过程为：[t, y] = ode45(myfun, [x0, xf], y0);其中myfun表示微分方程的函数形式，[x0, xf]表示求解的自变量范围，y0表示初值条件，t和y分别为求解得到的自变量和特解。
四、示例下面通过一个具体的示例来演示如何使用MATLAB求解微分方程特解。
假设我们要求解的微分方程为dy/dx = x^2 + y，初值条件为y(0) = 1，求解范围为x从0到5。
在MATLAB中定义微分方程的函数形式为：function dydx = myfun(x, y)dydx = x^2 + y;然后定义初值条件为：x0 = 0;y0 = 1;最后调用ode45函数求解微分方程特解：[t, y] = ode45(myfun, [x0, 5], y0);求解得到的自变量和特解分别存储在t和y中，可以通过绘图或其他方式对特解进行进一步分析。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下： function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵，D-L=tril(A,0) U=-triu(A,1); %求A的上三角阵 G=(D-L)\U; f=(D-L)\b; y=G*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end
第四讲 MATLAB解方程及微分方程 4.1 线性方程组求解 4.2 非线性方程数值求解 4.3 常微分方程初值问题的数值解法 4.4 函数极值
4.1 线性方程组求解
4.1.1 直接解法 1．利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运 算符“\”求解： x=A\b
例1 用直接解法求解下列线性方程组。
例2 用LU分解第2种格式求解例1中的线性方 程组。 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; %采用LU分解的第2种格式 [L,U ,P]=lu(A); x=U\(L\P*b)
Jacobi迭代法的MATLAB函数文件Jacobi.m如下： function [y,n]=jacobi(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵，注意diag(A)与diag(diag(A))区别 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵，注意-1与1的用处 B=D\(L+U); f=D\b; y=B*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=B*x0+f; n=n+1; end
例2 用LU分解求解例1中的线性方程组。 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [L,U]=lu(A); x=U\(L\b)
L= 1.0000 0 0 0 0.5000 1.0000 0 0 0 -0.3636 0.4773 1.0000 0.5000 -1.0000 1.0000 0 U= 2.0000 1.0000 -5.0000 1.0000 0 -5.5000 2.5000 6.5000 0 0 4.0000 2.0000 0 0 0 0.4091
例6 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初 值为0，迭代精度为10-6。
在命令中调用函数文件gauseidel.m，命令如下： A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]'; [x,n]=gauseidel(A,b,[0,0,0]',1.0e-6) x= 0.9958 0.9579 0.7916 n= 7
例4 用Cholesky分解求解例7中的线性方程组。 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; R=chol(A) ??? Error using ==> chol Matrix must be positive definite 命令执行时，出现错误信息，说明A为非正定矩阵。 讨论：[R,p]=chol(A)
例3 用QR分解第2种格式求解例1中的线性方 程组。 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; %采用QR分解的第2种格式 [Q,R,E]=qr(A); x=E*(R\(Q\b))
命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; %b=[13;-9;6;0]; x=A\b %x=inv(A)*b
x= -66.5556 25.6667 -18.7778 26.5556
2．利用矩阵的分解求解线性方程组
n= 4 y= NaN NaN NaN
m=
1012
2 非线性方程数值求解 2.1 单变量非线性方程求解
在MATLAB中提供了一个fzero函数，可以用 来求单变量非线性方程的根。该函数的调 用格式为： z=fzero('fname',x0,tol,trace) 其中fname是待求根的函数文件名，x0为搜索 的起点。一个函数可能有多个根，但fzero 函数只给出离x0最近的那个根。tol控制结 果的相对精度，缺省时取tol=eps，trace• 指 定迭代信息是否在运算中显示，为1时显示， 为0时不显示，缺省时取trace=0。
(2) QR分解
对矩阵X进行QR分解，就是把X分解为一个 正交矩阵Q和一个上三角矩阵R的乘积形式 (施密特正交化方法)。QR分解只能对方阵 进行。MATLAB的函数qr可用于对矩阵进 行QR分解，其调用格式为： [Q,R]=qr(X)：产生一个一个正交矩阵Q和一 个上三角矩阵R，使之满足X=QR。 [Q,R,E]=qr(X)：产生一个一个正交矩阵Q、 一个上三角矩阵R以及一个置换矩阵E，使 之满足XE=QR。 实现QR分解后，线性方程组Ax=b的解 x=R\(Q\b)或x=E(R\(Q\b))。
例3 用QR分解求解例1中的线性方程组。 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [Q,R]=qr(A); x=R\(Q\b)
Q= -0.8165 0 0.5173 -0.2564 -0.4082 -0.6848 -0.5955 0.0986 0 0.2490 -0.4301 -0.8678 -0.4082 0.6848 -0.4391 0.4142 R= -2.4495 -1.2247 4.4907 -2.0412 0 8.0312 -0.4358 -7.7822 0 0 -2.5775 -1.4647 0 0 0 -0.3550 Q'*Q= 1.0000 0 0.0000 -0.0000 0 1.0000 0.0000 -0.0000 0.0000 0.0000 1.0000 -0.0000 -0.0000 -0.0000 -0.0000 1.0000
2．Gauss-Serdel迭代法
在Jacobi迭代过程中，计算时，已经得到，不必再 用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改 进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比， Gauss-Serdel迭代用新分量代替旧分量，精度会 高些。
L= 1.0000 0 0 0 0.5000 1.0000 0 0 0.5000 -1.0000 1.0000 0 0 -0.3636 0.4773 1.0000 U= 2.0000 1.0000 -5.0000 1.0000 0 -5.5000 2.5000 6.5000 0 0 4.0000 2.0000 0 0 0 0.4091 P= 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0
Q= -0.1222 0.9576 0.0338 -0.2588 -0.8552 -0.0403 -0.4797 0.1920 0.1222 -0.1869 -0.5297 -0.8183 0.4887 0.2157 -0.6987 0.4759 R= -8.1854 0.2443 7.3302 -0.6108 0 -5.1904 2.0790 2.0905 0 0 -2.8189 -1.1108 0 0 0 0.1503 E= 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0
例5 用Jacobi迭代法求解下列线性方程组。设迭代初值为0， 迭代精度为10-6。
在命令中调用函数文件Jacobi.m，命令如下： A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]'; [x,n]=jacobi(A,b,[0,0,0]',1.0e-6) x= 0.9958 0.9579 0.7916 n= 11
(3) Cholesky分解 如果矩阵X是对称正定的，则Cholesky分解将矩阵
X分解成一个下三角矩阵和上三角矩阵的乘积。 设上三角矩阵为R，则下三角矩阵为其转置，即 X=R'R。MATLAB函数chol(X)用于对矩阵X进 行Cholesky分解，其调用格式为： R=chol(X)：产生一个上三角阵R，使R'R=X。若X 为非对称正定，则输出一个出错信息。 [R,p]=chol(X)：这个命令格式将不输出出错信息。 当X为对称正定的，则p=0，R与上述格式得到的 结果相同；否则p为一个正整数。如果X为满秩 矩阵，则R为一个阶数为q=p-1的上三角阵，且 满足R'R=X(1:q,1:q)。 实现Cholesky分解后，线性方程组Ax=b变成 R‘Rx=b，所以x=R\(R’\b)。
LU 分解又称Gauss 消去分解，可把任意方 阵分解为下三角矩阵的基本变换形式（行 交换）和上三角矩阵的乘积。即A=LU，L 为下三角阵，U 为上三角阵。线性代数中 已经证明，只要方阵A是非奇异的，LU分 解总是可以进行的。