代数方程和微分方程的计算机求解
微分方程计算机代数系中的应用
江西师范大学硕士学位论文计算机代数系统在微分方程研究中的应用姓名：曾广洪申请学位级别：硕士专业：基础数学指导教师：刘华祥20050401盐兰垫盛塾壅垫壅堕茎查堡堡墨±箜墨旦，ｇｉ：ｔ鞴ｐ工ｏｔ３ｄ‘【ｅｑｎ】ｔｆｘｔｔ），ｙ（ｔ’，ｚ（ｔ，】，ｔｔ０．．６００ｒ¨ｘ（Ｏ，互．２，ｙ《Ｏ）＃．２，ｚｆＯ，＝－３】ｊ，￥ｔｅｐｓｉｚｅ＝．０５，ｌｉｎｅｃｏｌｏｕｒ＝ｓｌｎ（ｔ＋ｐｉ／２））：＞９２：－ｔｔｘｔＰｌｏ七３ｄ（【１／（ａｌｐｈａ＋２），１／（ａｌｐｈａ＋２），Ｉ／（ａｌｐｈａ＋２），”一Ｐ”】，ｃｏｌｏｒ－ｒｅｄ，ｆｏⅡｔ－ＩＣＯＵＲＩＥＲ，ＢＯＬＤ，３０ｌ’：＞ｄｉｓｐｌａｙ（ｇｌ，９２，ｏｒｉｅｎｔａｔｉｏｎ＝【－２７，６０】，ｖｉｅｗ＝ｆｏ．２．．ｏ＋５，ｏ．２．．ｏ．５，ｏ．２．－Ｏ．５１）；圉６＿１更进～步，还可用以下语句给出动画图形，根据曲线随时间变化的趋势易验证系统在奇点Ｐ的邻域内有不稳定的周期解（因纸张不能呈现动画，故具体图形略去、．＞ｓｏｌ：。
ｄ８０ｌｖｅ《｛ｅｑｎ，ｘ（０）＊・２，了‘Ｏ）＊－２，２（Ｏ）－・３｝，｛ｘ‘ｔ），，ｆｔ），２（ｔ）＇，ｎｕｍｅｒｉｃ）：＞９３：－ｏｄｅｐｌ０２‘ｓｏｌ，ｆｘ（ｔ），ｙ（ｔ），＝《ｔ）ｌ…０６００，ｎｕｍｐｏｉａｔｓ－２４００，ｆｒａｍｅｓ－２０，ａｘｅｓｗＮＯＲＭＡＬ）：＞ｄｉｓｐｌａｙ（９２，９３ｒｏｒｉｅｎｔａｔｉｏｎ。
卜２７，６０】，ｖｉｅｗ＝【Ｏ．２．．０．５，０．２．．０．５，０．２．．０．ｓ¨；２４计算机代数系统在微分方程研究中的应用作者：曾广洪学位授予单位：江西师范大学本文链接：/Thesis_Y712281.aspx。
控制系统仿真_薛定宇第三章 科学运算问题的MATLAB求解
32/62
微分方程求解的步骤
将微分方程变换成标准型 用MATLAB描述微分方程
M-函数 入口：function dx=funmane(t,x) 匿名函数 >> f=@(t,x)[...]
求解
验证：odeset()函数
控制系统仿真与CAD 国家级精品课程
2014-12-31
2014-12-31
9/62
数值解
ห้องสมุดไป่ตู้
解析解
解析解
控制系统仿真与CAD 国家级精品课程
2014-12-31
10/62
演示：自编 funm() 求矩阵的任意函数
求
结果：左上角元素
控制系统仿真与CAD 国家级精品课程
2014-12-31
11/62
3.1 线性代数问题求解小结
线性代数很多问题可以用MATLAB语句直 接求解，和数学表示差不多一样直观 很多方法可以同时得出解析解和数值解
本节主要介绍和这门课程相关的问题 线性代数问题的MATLAB求解 代数方程求解、微分方程求解 最优化问题的求解 Laplace 变换与 z 变换问题的求解
控制系统仿真与CAD 国家级精品课程
2014-12-31
2/62
进一步学习这方面内容建议阅读
薛定宇、陈阳泉《高等应用数学问题的 MATLAB求解》（第二版），清华大学出 版社，2008。英文版：Solving Applied Mathematical Problems with MATLAB，CRC Press，2008 遇到某个MATLAB问题找不到合适的工具 箱，试在下面网址搜索
matlab求解最简单的一阶偏微分方程
matlab求解最简单的一阶偏微分方程一、引言在科学和工程领域，偏微分方程是非常重要的数学工具，用于描述各种现象和过程。
而MATLAB作为一种强大的数值计算软件，可以用来求解各种复杂的偏微分方程。
本文将以MATLAB求解最简单的一阶偏微分方程为主题，探讨其基本原理、数值求解方法以及具体实现过程。
二、一阶偏微分方程的基本原理一阶偏微分方程是指只含有一个未知函数的偏导数的微分方程。
最简单的一阶偏微分方程可以写成如下形式：\[ \frac{\partial u}{\partial t} = F(x, t, u, \frac{\partial u}{\partial x}) \]其中，\(u(x, t)\) 是未知函数，\(F(x, t, u, \frac{\partial u}{\partial x})\) 是给定的函数。
一阶偏微分方程可以描述很多实际问题，比如热传导、扩散等。
在MATLAB中，我们可以使用数值方法求解这类方程。
三、数值求解方法1. 有限差分法有限差分法是一种常用的数值求解偏微分方程的方法。
其基本思想是用离散的方式来逼近偏导数，然后将偏微分方程转化为代数方程组。
在MATLAB中，我们可以使用内置的求解器来求解离散化后的代数方程组。
2. 特征线法特征线法是另一种常用的数值求解方法，它利用特征线方程的特点来求解偏微分方程。
这种方法在求解一维情况下的偏微分方程时特别有效，可以提高求解的效率和精度。
四、MATLAB求解过程在MATLAB中，我们可以使用`pdepe`函数来求解一阶偏微分方程。
该函数可以针对特定的方程和边界条件，利用有限差分法进行离散化求解。
下面给出一个具体的例子来说明如何使用MATLAB求解最简单的一阶偏微分方程。
假设我们要求解如下的一维热传导方程：\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]其中，\(\alpha\) 是热传导系数。
线性代数与计算机的关系
线性代数与计算机的关系线性代数与计算机的关系————————————————————————————————作者：————————————————————————————————日期：浅谈线性代数与计算机的关系线性代数是计算机专业的一门重要基础课程,同时又作为各高等院校和工科类专业的数学基础课程，它具有很强大的应用性和实用性。
线性代数是数学的一个分支，它主要处理线性关系问题，它的研究对象是向量、向量空间、线性变换和有限维的线性方程组，向量空间是现代数学的一个重要课题；因而，线性代数被广泛应用于抽象代数和泛函分析中；用过解析几何，线性代数得以被具体表示。
线性代数的理论已经被泛化为算子理论。
由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。
自计算机产生以来，随着计算机的不断发展和进步，计算机语言也在进步，但是很多软件或编程的编写都离不开计算机算法，这时一种好的计算方法就会成为一个软件或编程的亮点。
以前，在计算机的计算算法中，对于一些复杂的计算总是要花很多步骤来完成，既麻烦又容易出错，并很浪费时间（比如在计算机上用算法求鸡兔同笼的问题，如果是用一般算法来求的话，我们会发现很吃力，但是引用的线性代数的矩阵理论就简单的多了），所以在计算效率方面提不上去的话，就会限制计算机的发展和进步。
而线性代数的引入就改变了这个问题，使得计算机的发展更加迅猛，到了今天计算机得到广泛应用的时候，计算机数据结构、算法、计算机图形学、计算机辅助设计、密码学、经济学、网络技术、虚拟现实等技术无不是以线性代数为理论基础并组成其计算机算法中极其重要的一部分。
线性代数在计算机领域的应用与计算机的计算性能是成正比例的，同时，这一性能会随着计算机硬件的不断创新和发展而得到极大的提升。
线性代数的计算机应用在全球有很多的应用，例如Wassily Leontief教授把美国经济用500个变量的500个线性方程组描述,而后又把系统简化为42个变量的42个线性方程。
第8章 代数方程和常微分方程求解
8.2 常微分方程求解
求解微分方程必须事先对自变量的某些值规定出 函数或是导数的值。 若在自变量为零的点上，给出初始条件，称为初 值问题，最普遍的自变量是“时间”。例如，弹 性系统的自由振动，若以时间为零来限定位移和 速度，这是一个初值问题。 若在自变量为非零的点上，给出边界条件，称为 边值问题，最普遍的自变量是“位移”。例如， 描述梁弯曲变形的微分方程，边界条件总是规定 在梁的两端。
当 x 0 2 和 y 0 0 条件下的特解。 在此问题中，两个微分方程的MATLAB表达式为： e1：Dx+2*x-Dy=10*cos(t) e2：Dx+Dy+2*y=4*exp(-2*t) 初值条件表达式为： C1：x(0)=2 C2：y(0)=0
8.1 代数方程求解
8.1.1 代数方程图解法
符号绘图函数fplot()和ezplot()也可以用于图解 法求代数方程的根，它适用于求解维数较少的一 维方程或二维方程组。 对于一维方程图解，其解就是函数曲线与x轴交点 所对应的变量数值。如果有多个交点，则表示该 方程有多个解；如果没有交点，则表示该方程没 有解。 例如，在例5-3使用符号绘图函数绘制代数方程的 图形（图5-3左图）中可见，函数在区间[-5,5]内 与x轴有3个交点，因此该代数方程该区间内有3个 实根。
M文件运行结果： 采用矩阵左除或矩阵求逆求出线性方程组的解： xx (zx)= 1.0000 2.0000 3.0000 -1.0000 计算残量： r = 1.0e-014 * 0.0888 0.2220 -0.4441 0.1776 计算残量的模： R = 5.3475e-015
偏微分方程的解法
偏微分方程的解法偏微分方程（Partial Differential Equations，简称PDEs）是数学中的一个重要分支，它描述了多变量函数的偏导数之间的关系。
这些方程在自然科学、工程应用和社会科学等领域都发挥着重要作用。
解决偏微分方程是一个复杂而有挑战性的过程，需要运用多种数学方法和工具来求解。
在本文中，我将为您介绍几种常见的偏微分方程的解法，并提供一些示例以帮助您更好地理解。
以下是本文的主要内容：1. 一阶线性偏微分方程的解法1.1 分离变量法1.2 特征线方法2. 二阶线性偏微分方程的解法2.1 分离变量法2.2 特征值法2.3 Green函数法3. 非线性偏微分方程的解法3.1 平移法3.2 线性叠加法3.3 变换法4. 数值方法解偏微分方程4.1 有限差分法4.2 有限元法4.3 谱方法5. 偏微分方程的应用领域5.1 热传导方程5.2 波动方程5.3 扩散方程在解一阶线性偏微分方程时，我们可以使用分离变量法或特征线方法。
分离变量法的基本思路是将方程中的变量分离，然后通过积分的方式求解每个分离后的常微分方程，最后再将结果合并。
特征线方法则是将方程中的变量替换为新的变量，使得方程中的导数项消失，从而简化求解过程。
对于二阶线性偏微分方程，分离变量法、特征值法和Green函数法是常用的解法。
分离变量法的核心思想与一阶线性偏微分方程相似，将方程中的变量分离并得到常微分方程，然后进行求解。
特征值法则利用特征值和特征函数的性质来求解方程，适用于带有齐次边界条件的问题。
Green函数法则通过引入Green函数来求解方程，其特点是适用于非齐次边界条件的情况。
非线性偏微分方程的解法则更加复杂，常用的方法有平移法、线性叠加法和变换法。
这些方法需要根据具体问题的特点选择合适的变换和求解技巧，具有一定的灵活性和创造性。
除了上述解析解法，数值方法也是解偏微分方程的重要手段。
常用的数值方法包括有限差分法、有限元法和谱方法等。
微分方程的经典解法
01
02
03
非线性变量代换法
变量代换法的应用
变量代换法在解决各种实际问题中有着广泛的应用，如物理、工程、经济等领域。
通过选择适当的代换变量，可以简化复杂的微分方程，从而更方便地求解。
变量代换法是解决微分方程的一种重要技巧，尤其在处理非标准形式的微分方程时非常有效。
01
高阶非线性微分方程的解法通常包括迭代法、摄动法和数值方法等。
02
迭代法是通过不断迭代方程的解来逼近真实解，常用的方法有牛顿迭代法和欧拉迭代法等。
03
摄动法是将非线性微分方程转化为摄动方程，然后通过小参数展开求解。
04
数值方法是通过离散化微分方程，然后使用计算机求解离散化后的方程组。
高阶微分方程在物理、工程、经济等领域有广泛应用，如振动分析、控制系统、信号处理等。
04
积分因子法
积分因子法是一种求解微分方程的方法，通过引入一个积分因子来消除方程中的导数项，从而将微分方程转化为代数方程进行求解。
积分因子法适用于可分离变量、线性、部分线性以及某些非线性微分方程。
积分因子法的关键是找到一个函数，使得该函数与微分方程的每一项相乘后，能够消去方程中的导数项。
方法概述
高阶线性微分方程的一般形式为$y^{(n)}(x) + a_{n-1}(x)y^{(n-1)}(x) + cdots + a_0(x)y(x) = 0$。
变量分离法是将方程转化为多个一阶微分方程，然后分别求解。
幂级数法是通过将解表示为幂级数的形式，然后代入初始条件求解系数。
高阶非线性微分方程的解法
02
通过引入新变量 (u = ax + by)，可以将原方程转化为 (y^{prime} = frac{1}{a} f(u))。
微分方程式的建立与求解
通过建立微分方程式描述物体在重力作用下的运动规律，如速度、加速度与时 间的关系。
02
微分方程的求解方法
分离变量法
总结词
通过将微分方程转化为代数方程，简 化求解过程。
详细描述
分离变量法适用于具有两个变量的微 分方程，通过分离变量，将微分方程 转化为代数方程，然后求解代数方程 得到微分方程的解。
05
微分方程的稳定性分析
线性微分方程的稳定性分析
线性微分方程的稳定性分析主要基于其 特征值和特征向量。如果所有特征值都 位于复平面的左半部分，则系统是稳定 的；否则，系统是不稳定的。
线性微分方程的解可以通过求解其特征值和 特征向量得到，也可以通过积分得到。
线性微分方程的解具有叠加性，即 如果两个解都是稳定的，那么它们 的线性组合也是稳定的。
振动分析
在研究物体的振动时，通过建立位移、速度和加 速度的微分方程来分析振动的规律和特性。
3
热传导方程
在研究热量在物体中的传递时，通过建立温度关 于时间和空间的微分方程来模拟热传导过程。
在经济中的应用
供需关系
01
在分析商品市场的供需关系时，通过建立需求和供给函数的微
分方程来预测价格变动。
经济增长模型
非线性微分方程的稳定性分析
非线性微分方程的稳定性分析比线性微分方程更为复杂，需要考虑更多的因素，如非线性项的性质、 初始条件等。
非线性微分方程的解可以通过数值方法（如欧拉法、龙格-库塔法等）得到，也可以通过解析方法（如 分离变量法、幂级数展开等）得到。
非线性微分方程的解具有不可叠加性，即如果两个解都是稳定的，那么它们的线性组合不一定是稳定的。
微分方程式的建立与 求解
目 录
《数学实验》课程教学大纲
3．线性方程组 的求解
4.随机数的生产和模拟
5．实验实例：循环比赛的名次和按年龄分组的种最优化问题实验
重点：学会一些常用函数的调用格式并学会自己动手编写函数
3. 《高等应用数学问题的MATLAB求解》.薛定宇,陈阳泉著.清华大学出版社,2004
4. 《MATLAB数学实验》.胡良剑,孙晓君编著.高等教育出版社，2006.6
执笔人：邓化宇
审核人：
院（系）负责人：
《数学实验》课程教学大纲
MathematicalExperiment
适用：本科四年制信息与计算科学专业(40学时左右)
一、课程的目的及任务
开设《数学实验》课的目的是在两周的时间里为学生介绍如何使用计算机的语言和方法去处理一些经典的数学问题，并提供一些实例以启发学生自己动手练习。进一步的提高要靠学生的兴趣和努力。
教学要点：
1．一元非线性方程数值求解
2．非线性方程组数值求解
3．方程符号求解
4．一元函数和多元函数无约束优化求解
5．线性规划
6．实验实例：购房贷款的利率和最短路问题
第五章 微分方程问题的计算机求解
重点：学会一些常用函数的调用格式并学会自己动手编写函数
教学要点：
1．常系数微分方程的计算机求解析解
2．微分方程问题的数值解法
二、课程的特点、要求及本课程与其它课程的联系
数学是科学技术人才科学素质的的重要组成部分，随着高科技与与计算技术的发展和普及，数学的重要性日益突出。“高技术本质上是一种数学技术”这一观点已越来越多地为人们所认同。学习计算机使用和开发是启迪学生创新意识和创新思维、锻炼创新能力、培养高层次人才的一条重要途径；也是激发学习欲望、培养主动探索、努力进取学风和团结协作精神的有力措施。
实验报告—代数方程与微分方程求解
实 验 报 告 四代数方程求解1、【示例】以下命令可求出方程 (x +1)e –x +e x sin x =0在0附近的一个根：>>y=sym('(x+1)*exp(-x)+exp(x)*sin(x)'); % 用sym 命令定义符号表达式>>x=solve(y,'x') % 用准解析方法求出方程最接近0的一个根x =-0.86508244315736795185621568221837或可用以下命令求解该方程以指定点为初始搜索点的数值解：>> y=inline('(x+1)*exp(-x)+exp(x)*sin(x) ', 'x'); % 用数值方法求解时，方程要用inline 命令定义 >> x=fsolve(y,0) % 用数值方法从初始点1开始搜索方程的近似解 x = -0.8651注：准解析命令solve 只能求出方程最接近0的一个实数根，而数值解法fsolve 可以通过初始搜索点的变化，得到不同的解（如果方程有多个实数解）。
【要求】仿照示例，用准解析方法求出30.5sin(42)4cos(2)0.5t t e t e t --++=的一个根；再用数值方法分别求该方程在-0.6和3附近的两个根。
y=sym('exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5');t=solve(y,'t')t =0.67374570500134756702960220427474y=inline('exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5','t');t=fsolve(y,0.6)t =0.6737y=inline('exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5','t');t=fsolve(y,3)t =2.59372、【示例】以下命令可求解非线性方程组339820x y x x y ⎧+-=⎨+-=⎩>> eq1=sym('x^3+y^3-x-98'); % 定义第一个方程表达式>> eq2=sym('x+y-2'); % 定义第二个方程表达式>> [x,y]=solve(eq1,eq2) % 解方程组（用准解析方法）x =13/12+1/12*2329^(1/2)13/12-1/12*2329^(1/2)y =11/12-1/12*2329^(1/2)11/12+1/12*2329^(1/2)或可用以下命令求解上述方程组以指定点为初始搜索点的数值解：>> f=inline('[x(1) ^3+x(2) ^3-x(1)-98; x(1)+x(2)-2]', 'x'); % 用inline 命令定义方程组>> x=fsolve(f,[1;1]) % 用数值方法从初始点(1,1)开始搜索方程组的一个近似解 x =-2.93834.9383【要求】仿照示例，求解（1）方程35323=+-x x x 和方程组⎪⎩⎪⎨⎧-==+=++1430122yz z x x x 35323=+-x x x（1）y=sym('x^3-3.*x^2+5*x');x=solve(y,'x')x =1.1.+1.4142135623730950488016887242097*i1.-1.4142135623730950488016887242097*i（2）y=inline('x^3-3.*x^2+5*x-3 ', 'x');x=fsolve(y,0)x =1.0000⎪⎩⎪⎨⎧-==+=++1430122yz z x x x（1）eq1=sym('x^2+2*x+1');eq2=sym('x+3*z-4');eq3=sym('y*z+1');[x,y,z]=solve(eq1,eq2,eq3)x =-1y =-3/5z =5/3（2）f=inline('[x(1)^2+2*x(1)+1; x(1)+3*x(3)-4;x(2)*x(3)+1]', 'x');x=fsolve(f,[1;1;1])x = -0.9979-0.60031.6660（2）方程组x 2+y =0, e x +sin y =0在[-5,5]内的所有根。
数值分析基础
数值分析基础数值分析是一门研究利用计算机进行数值计算的学科，它涉及到数学、计算机科学和工程学等多个领域。
数值分析基础是数值计算领域最基本的理论和方法，为实现高精度、高效率的数值计算提供了重要的基础。
一、数值分析的概念数值分析是通过数值方法解决数学问题的过程。
它的基本思想是将连续的数学问题转化为离散的数值问题，并利用计算机进行求解。
数值分析的应用范围非常广泛，包括线性代数方程组的求解、非线性方程求根、插值与逼近、数值微积分、常微分方程的初值问题和边值问题的数值解等。
二、数值计算的误差分析在数值分析中，误差分析是非常重要的一环。
数值计算过程中产生的误差可以分为截断误差和舍入误差。
截断误差是由于在离散化和近似计算中引入的近似误差，而舍入误差是由于计算机在表示实数时的有限精度引起的。
准确估计和控制误差是数值计算的核心问题之一。
三、常用的数值计算方法1. 插值与逼近方法：插值是在给定一组数据点的情况下，通过构造一个函数来近似这组数据点之间未知函数值的方法。
常用的插值方法有拉格朗日插值和牛顿插值。
逼近是通过在给定函数空间中寻找一个尽可能接近原函数的近似函数的方法，常见的逼近方法有最小二乘逼近和Chebyshev逼近。
2. 数值积分方法：数值积分是计算定积分的近似值的方法。
常见的数值积分方法有梯形法则、辛普森法则和复合求积法。
3. 数值微分方法：数值微分是通过差商逼近导数的计算方法。
常见的数值微分方法有中心差商、前向差商和后向差商。
4. 数值求解线性方程组的方法：线性方程组求解是数值计算中的一个重要问题。
常用的求解方法有直接法和迭代法。
5. 常微分方程数值解法：常微分方程数值解法是通过数值方法求解微分方程的方法。
常用的数值解法有欧拉法、龙格-库塔法和变步长方法等。
四、数值计算的应用领域数值分析在各个学科领域都有广泛的应用。
在物理学中，数值分析被用于求解天体运动、弹道问题等。
在工程学中，数值分析被用于优化设计、结构力学分析等。
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)单击工具栏的“=”按钮，开始求解。
微分代数方程
微分代数方程(DAE)的Matlab解法所谓微分代数方程，是指在微分方程中，某些变量满足某些代数方程的约束。
假设微分方程的更一般形式可以写成前面所介绍的ODEs数值解法主要针对能够转换为一阶常微分方程组的类型，故DAE就无法使用前面介绍的常微分方程解法直接求解，必须借助DAE的特殊解法。
其实对于我们使用Matlab求解DAE时，却没有太大的改变只需增加一个Mass参数即可。
描述f(t,x)的方法和普通微分方程完全一致。
注意：ode15i没法设置Mass属性，换句话说除了ode15i外其他ode计算器都可以求解DAEs问题1.当M(t,y)非奇异的时候，我们可以将微分方程等效的转换为y'=inv(M)*f(t,y)，此时就是一个普通的ODE(当然我们可以将它当成DAEs处理)，对任意一个给定的初值条件都有唯一的解2.当m(t,y)奇异时，我们叫它为DAEs(微分代数方程)，DAEs问题只有在同时提供状态变量初值y0和状态变量一阶导数初值py0，且满足M(t0,y0)*yp0=f(t0,y0)时才有唯一解，假如不满足上面的方程，DAEs解算器会将提供的y0和py0作为猜测初始值，并重新计算与提供初值最近的封闭初值3.质量矩阵可是一个常数矩阵(稀疏矩阵)，也可以是一个自定义函数的输出。
但是ode23s只能求解Mass是常数的DAEs4.对于Mass奇异的DAEs问题，特别是设置MassSingular为yes时，只能ode15s 和ode23t解算器5.对于DAE我们还有几个参数需要介绍a.Mass：质量矩阵，不说了，这个是DAE的关键，后面看例子就明白了b.MStateDependence：质量矩阵M(t,y)是否是y的函数，可以选择none|{weak}|strong，none表示M与y无关，weak和strong都表示与y相关c.MvPattern：注意这个必须是稀疏矩阵，S(i,j)=1表示M(t,y)的第i行中任意元素都与第j个状态变量yi有关，否则为0d.MassSingular：设置Mass矩阵是否奇异，当设置为yes时，只能使用ode15s 和ode23te.InitialSlope：状态变量的一阶导数初值yp0，和y0具有相同的size，当使用ode15s和ode23t时，该属性默认为0下面我们以实例说明，看下面的例子，求解该方程的数值解【解】真是万幸，选取状态变量和求状态变量的一阶导数等，微分方程转换工作，题目已经帮我们完成。
代数方程和微分方程求解PPT教学课件
2020/12/12
8
多项式运算的几个常用函数：
P=conv(p1,p2);
%多项式乘法
[d,r]=deconv(p1,p2); %多项式除法
Dp=polyder(p);
%多项式的导数
Ip=polyint(p) %多项式的积分（原函数）
Y=polyval(p,x)
%输出多项式p在向量x
的值
2020/12/12
9
多项式求根的函数为 r=roots(p)
求得多项式p的所有根。
例4.3：求多项式 x 6 2 x 5 0 1 x 4 3 3 x 8 3 2 2 x 8 2 2 13 x 6 19 22 6 的根并在多项式图形中表示。
2020/12/12
10
参考程序：
q =[1 -20 138 -328 -223 1692 1260]; r=roots(q); x=-2.2:0.05:8; y=polyval(q,x); y1=polyval(q,r); plot(x,y,r,y1,'p') xlim([-2.2,8]) legend('polynomial','roots')
2020/12/12
11
2020/12/12
12
线性方程组的求解
线性方程组 Ax=b
可以利用矩阵除法直接得到。但当系数矩阵为稀 疏矩阵时，利用稀疏矩阵函数可以得到更高的计 算效率。 稀疏矩阵利用函数
A1=sparse(A); 定义。
2020/12/12
13
例4.4：求解n阶线性方程组
2 1
0 x1 1
4
也可以利用下面的语句求解
>>
基于Matlab常系数线性微分方程组的求解
基于Matlab常系数线性微分方程组的求解严水仙【摘要】在常微分方程课程教学中,常系数线性微分方程组可以通过线性代数的理论、矩阵指数、拉普拉斯变换等方法进行求解.本文主要叙述利用Matlab数学软件在求解常系数线性微分方程组中的应用.【期刊名称】《赣南师范学院学报》【年(卷),期】2018(039)003【总页数】5页(P10-14)【关键词】常系数线性微分方程;Matlab;矩阵指数【作者】严水仙【作者单位】赣南师范大学数学与计算机科学学院,江西赣州341000【正文语种】中文【中图分类】O175微分方程课程是高校不少理工科专业(如数学、力学、控制等) 的重要基础理论课程.常微分方程是描述自然科学、工程技术和社会科学中的运动、演化和变化规律的重要连续型模型. 物理、化学、材料、医学、经济学等领域中的许多原理和规律都可以描述成相应的微分方程, 如生物种群中的生态平衡、流行病存在的阈值定理、化学反应中的稳定性、遗传基因变异、股票的涨幅趋势、利率的浮动、市场均衡价格的变化等.描述、认识和分析其中的规律可以通过研究相应的微分方程数学模型来实现.[1]在微分方程的理论中，线性微分方程组是非常值得重视的一部分内容，它是了解并掌握非线性微分方程、非线性动力系统、非线性控制等课程的基础. 常系数线性微分方程组的求解是线性微分方程组理论中最简单、最直观的部分，熟悉并掌握常系数线性微分方程的求解将有利于更好的理解线性系统的基本理论.Matlab是由美国的Cleve Moler博士等[2-3]于1980年提出的以矩阵运算为基础，把计算、程序设计等融合到了一个简单易用的交互式工作环境中.可实现工程计算、算法研究、符号运算、建模和仿真、原型开发、数据分析及可视化、科学和工程绘图、应用程序设计等功能.Matlab强大的运算功能和图形使其成为目前世界上应用最为广泛的科学计算软件之一, 在教学中能快速的计算方程的解并描绘直观的几何图形.[4-6]鉴于此,本文主要介绍借助于Matlab来求解常系数线性微分方程组，通过利用Matlab命令，计算系数矩阵的特征值、特征向量、矩阵指数求解线性微分方程组.1 常系数线性微分方程的基本理论[1]定理1[1] 如果A(t)是n×n阶矩阵函数，f(t)是n维列向量函数.它们都在区间a≤t≤b上连续，则对区间a≤t≤b上的任意t0∈[a,b]及任一常数n维列向量η,方程组x′=A(t)x+f(t)(1)存在唯一解φ(t)，定义于整个区间a≤t≤b上，且满足初值条件φ(t0)=η.定理2[1] 齐次线性微分方程组x′=A(t)x一定存在n个线性无关的解x1(t),x2(t),…,xn(t).定理 3[1] 齐次线性微分方程组x′=A(t)x一定存在一个基解矩阵Φ(t).如果ψ(t)是方程组的任意解，那么ψ(t)=Φ(t)c,(2)这里c是确定的n维常数列向量.定理4[1] 如果Φ(t)，ψ(t)在区间a≤t≤b上是x′=A(t)x的两个基解矩阵，那么，存在一个非奇异n×n常数矩阵C，使得在a≤t≤b区间上ψ(t)=Φ(t)C.定理5[1] 设Φ(t)是齐次线性微分方程组x′=A(t)x的基解矩阵，是非齐次线性微分方程组x′=A(t)x+f(t)的某一个解，则方程组x′=A(t)x+f(t)的任一解φ(t)都可表示为(3)这里c是确定的n维常数列向量.矩阵指数exp A的定义和性质：如果A是一个n×n常数矩阵，则定义为矩阵指数，其中E为n阶单位矩阵，且规定A0=E,0!=1,对于所有元均为0的矩阵，易知expO=E.根据矩阵指数的定义及级数的收敛性，易知对一切矩阵A都是绝对收敛，在t的任何有限区间上是一致收敛.矩阵指数expA有如下性质：(Ⅰ)如果矩阵A,B是可交换的矩阵，即AB=BA,则exp(A+B)=exp A exp B. (Ⅱ)对于任何矩阵A,(exp A)-1存在，且(exp A)-1=exp(-A).(Ⅲ)如果T是非奇异矩阵(可逆矩阵)，则exp(T-1AT)=T-1(exp A)T.设常系数线性微分方程组为x′=Ax(4)其中A是n×n阶常数矩阵,x=(x1,x2,…,xn)T.定理6 矩阵Φ(t)=exp At(5)为方程组(4)的基解矩阵，且Φ(0)=E.证明由矩阵指数的定义易得Φ(0)=E.(1.5)式对t求导，我们得到Φ′=(exp A t) exp At=AΦ(t),所以Φ(t)=exp At是方程(4)的解矩阵.又因为det Φ(0)=det E=1,因此Φ(t)是(4)的基解矩阵.定理6已经给出了常系数线性微分方程组(4)的基解矩阵，但是exp At是由At的矩阵级数定义的，矩阵中的每个元是什么没有具体给出，下面我们讨论exp At的计算方法.2 常系数线性微分方程组基解矩阵的求解定理7 如果矩阵A具有n个线性无关的特征向量v1,v2,…,vn，它们对应的特征值分别为λ1,λ2,…,λn(不必各不相同)，那么矩阵ψ(t)=[eλ1tv1,eλ2tv2,…,eλntvn],(-∞≤t≤+∞)是常系数线性微分方程组(4)的一个基解矩阵.证明因为λ1,λ2,…,λn为矩阵A的特征值，对应的特征向量为v1,v2,…,vn，有det(λiE-A)=0,且(λiE-A)vi=0,i=1,2,…,n.又因为eλit≠0，则λieλitvi=Aeλitvi，即φi(t)=eλitvi是方程组(4)的一个解.因此，矩阵ψ(t)=[eλ1tv1,eλ2tv2,…,eλntvn]是(4)的一个解矩阵.v1,v2,…,vn是矩阵A线性无关的特征向量组，所以det ψ(0)=det[v1,v2,…,vn]≠0，因此可知，ψ(t)=[eλ1tv1,eλ2tv2,…,eλntvn]是(4)的一个基解矩阵.由此可知，Φ(t)=exp At和ψ(t)=[eλ1tv1,eλ2tv2,…,eλntvn]均为方程组(4)的基解矩阵，根据定理4的结论可知，存在一个非奇异的常数矩阵C,使得Φ(t)=expAt=ψ(t)C.令t=0，我们得到C=ψ-1(0)，因此exp At=ψ(t)ψ-1(0).所以，常系数线性微分方程的求解转化为求系数矩阵的特征值和特征向量.例1 求下列方程组的基解矩阵解系数矩阵为的特征方程为系数矩阵的特征值为λ1=1,λ2=2,λ3=3.设λ1=1时对应的特征向量v1=(a,b,c)T，则(λ1E-A)v1=0，解得特征向量为其中c1为参数.设λ2=2时对应的特征向量v2=(a,b,c)T，则(λ2E-A)v2=0，解得其中c2为参数. 设λ3=3时对应的特征向量v3=(a,b,c)T，则(λ3E-A)v3=0，解得其中c3为参数.由定理7知，方程组的基解矩阵为：根据定理3，方程组的通解为φ(t)=ψ(t)C,其中C=(c1,c2,c3)T.常系数线性微分方程的求解实际就是求系数矩阵的特征值和特征向量，或者通过拉普拉斯变换法求解以及直接求解系数矩阵的矩阵指数exp(At). 然而，不管矩阵的的特征值、特征向量、拉普拉斯变换、还是矩阵指数，直接求解均比较复杂，特别是系数矩阵的特征值出现重根、复根等情况时，求解通解就显得更为困难.下面介绍使用Matlab数学软件求解常系数线性方程组.3 Matlab在求解常系数微分方程中的应用根据线性代数、高等代数的理论，矩阵对角化过程是一个复杂的计算过程，特别是系数矩阵的特征值出现重根、复根等情况时，特征向量计算比较麻烦. 学习使用Matlab软件在求解数学问题的原理将能够让学生更好的理解数学思想，减少重复性计算等. 下面将介绍几个利用Matlab计算基解矩阵的例子.例2 求下列线性微分方程组满足初始条件φ(0)=η的解其中解在Matlab软件中直接输入如下命令，>> A=[0 1 0;0 0 1;-6 -11 -6];>> syms t;>> expm(A*t) %运行得线性方程组的基解矩阵exp(At)ans =[3*exp(-t) - 3*exp(-2*t) + exp(-3*t), (5*exp(-t))/2 - 4*exp(-2*t) + (3*exp(-3*t))/2, exp(-t)/2 - exp(-2*t) + exp(-3*t)/2][6*exp(-2*t) - 3*exp(-t) - 3*exp(-3*t), 8*exp(-2*t) - (5*exp(-t))/2 - (9*exp(-3*t))/2, 2*exp(-2*t) - exp(-t)/2 - (3*exp(-3*t))/2][3*exp(-t) - 12*exp(-2*t) + 9*exp(-3*t), (5*exp(-t))/2 - 16*exp(-2*t) +(27*exp(-3*t))/2, exp(-t)/2 - 4*exp(-2*t) + (9*exp(-3*t))/2]即得方程组的基解矩阵所以方程组的通解为φ(t)=exp(At)C,其中C=(c1,c2,c3)T.初始条件φ(0)=η代入φ(t)=exp(At)C，易得C=(1,2,2)T，故满足初始条件的解为例3 求下列线性微分方程组的通解解在Matlab软件中直接输入如下命令，>> A=[3 -1 1;2 0 1;1 -1 2];>> syms t;>> expm(A*t) %运行得线性方程组的基解矩阵exp(At)ans =[exp(2*t) + t*exp(2*t), - exp(2*t) - t*(exp(2*t) - exp(2*t)/t), exp(2*t) +t*(exp(2*t) - exp(2*t)/t)][exp(2*t) - exp(t) + t*exp(2*t), exp(t) - exp(2*t) - t*(exp(2*t) - exp(2*t)/t), exp(2*t) + t*(exp(2*t) - exp(2*t)/t)][exp(2*t) - exp(t), exp(t) - exp(2*t), exp(2*t)].即得方程组的基解矩阵所以方程组的通解为φ(t)=exp(At)C,其中C=(c1,c2,c3)T. 由此可知，通过使用Matlab软件求解常系数线性微分方程组的通解变得简单、直观.在常微分方程的理论和应用上，微分方程组的解在t→∞时解的性态的研究是很重要的内容，其对非线性方程的理论有着非常重要的应用.4 常系数线性方程组的稳定性定理定理8 常系数线性微分方程组(4)的解有下列性质:(Ⅰ)如果系数矩阵A的特征值的实部都是负的，则(4)的任一解当t→∞时都趋于零(稳定)；(Ⅱ)如果系数矩阵A的特征值的实部至少有一个是正的，则(4)至少有一个解t→∞时趋于无穷(不稳定)；(Ⅲ)如果系数矩阵A的特征值的实部都是非正，但有零根或具有零实部的根，则(4)的解在t→∞时不能判定(不确定).根据定理8，很容易判定常系数线性微分方程组解的性态.当t→∞时，例1、3的解趋于无穷，例2的解趋于零.由Matlab软件也可直接描述出方程组解的变化趋势.下面画出例2的解曲线图.在Matlab中直接输入命令：>> t=0∶200;>>plot(t,8*exp(-1*t)-13*exp(-2*t)+5*exp(-3*t),′o′,t,-8*exp(-1*t)+26*exp(-2*t)-15*exp(-3*t),′*′,t,8*exp(-1*t)-52*exp(-2*t)+45*exp(-3*t),′k′) %运行即得解在t∈[0,200]随t的变化图.5 小结常系数线性微分方程的求解是常微分方程理论中重要的组成部分，但计算比较复杂,对学生的数学基础要求较高,所以学生学起来比较吃力,容易对学习丧失兴趣.所以,在授课过程中,我们不仅要将基本概念和原理给学生讲通讲透,重点介绍数学思维，数学思想，还要利用计算机的表现能力将抽象的问题具体化，复杂的计算简单化.Matlab教学平台的引入,能够化繁为简,化抽象为具体,加深学生对本课程的掌握程度,提高教学效果,并且引导学生将理论应用于实际.【相关文献】[1] 王高雄,朱思铭,等.常微分方程[M].第3版.北京：高等教育出版社，2006.[2] 胡良剑,孙晓君.MATLAB数学实验[M].第2版.北京：高等教育出版社，2014.[3] 朱春蓉,郑群珍.Maple在常微分方程教学中的应用[J].河南教育学院学报：自然科学版，2009,18(3):63-64.[4] 闫金亮.Matlab在常微分方程教学中的应用[J].武夷学院学报，2012,31(2):96-100.[5] 张守贵.一类二阶常系数微分方程特解的教学探讨[J].重庆工商大学学报：自然科学版，2012, 29(12):12-15.[6] 徐定华,葛美宝.论微分方程课程的教学设计[J].大学数学, 2010,26(3):1-5.。
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)单击工具栏的“=”按钮，开始求解。
第12章 MATLAB常微分方程(组)数值求解方程与方程组的数值解
ode15s求解起来效率不太高的刚性问题。
ode23t可以用来求解微分代数方程。
ode23tb 刚性
低
当方程是刚性的，并且求解要求精度不高
时可以使用。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
常微分方程数值解
【例12.4-1】的结果图：
方 法 1计 算 结 果 图 1
方 法 2计 算 结 果 图 1
方 法 3计 算 结 果 图 1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0
0
0
-0.2
-0.4 0
2020/6/19
y1(t)
-0.2
y2(t)
y3(t)
-0.4
在某段时间区间内的变化来看。非刚性问题变化相对缓 慢，而刚性问题在某段时间内会发生剧烈变化，即很短 的时间内，解的变化巨大。对于刚性问题不适合用 ode45来求解，如果硬要用ode45来求解的话，达到指定 精度所耗费的时间往往会非常长 。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
二、 非刚性问题举例
这些函数可以求解非刚性问题，刚性问题，隐式
微分方程，微分代数方程等初值问题，也可以求解 延迟微分方程，以及边值问题等。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
二、初值问题求解函数
常微分方程数值解
1. 提供的函数
ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb，这些函数统一 的调用格式如下：
matlab有限元求解微分方程的本征值
一、概述Matlab是一种常用的数学软件，它提供了丰富的工具和函数，可用于解决各种数学问题。
其中，有限元法是一种常用的数值求解方法，它可用于求解微分方程的本征值问题。
本文将探讨如何使用Matlab进行有限元求解微分方程的本征值问题。
二、有限元法简介有限元法是一种数值分析方法，它通过将连续的物理问题离散化为有限数量的单元或网格，然后利用线性代数方法求解离散问题，从而得到原始的连续问题的近似解。
在微分方程的求解中，有限元法可用于求解微分方程的本征值问题，即确定微分方程的本征值和本征函数。
三、使用Matlab进行有限元求解微分方程的本征值问题1. 离散化微分方程需要将微分方程离散化为有限元形式。
这通常涉及将微分方程转化为一个矩阵形式的代数方程组。
对于一维问题，可以将区域离散化为一系列节点，并将微分方程表示为每个节点上的代数方程。
对于二维或三维问题，可以将区域离散化为网格或单元，并在每个单元中求解微分方程。
2. 构建刚度矩阵和质量矩阵一旦微分方程被离散化，就可以构建刚度矩阵和质量矩阵。
刚度矩阵描述了系统的刚度和连接性，质量矩阵描述了系统的质量和惯性。
这两个矩阵可以通过有限元方法和数值积分计算得到。
3. 求解本征值问题一旦刚度矩阵和质量矩阵被构建，就可以通过求解本征值问题来得到微分方程的本征值和本征函数。
这通常涉及求解特征值问题，即寻找一个非零向量，使得矩阵乘以该向量等于特征值乘以该向量。
4. 使用Matlab进行求解Matlab提供了丰富的工具和函数，可用于构建刚度矩阵和质量矩阵，并求解本征值问题。
使用Matlab的有限元工具箱或相关函数，可以方便地进行有限元求解微分方程的本征值问题。
四、案例分析下面通过一个简单的例子来说明如何使用Matlab进行有限元求解微分方程的本征值问题。
考虑一维弦的振动问题，其微分方程为：$$\frac{d^2u}{dx^2} +\omega^2u = 0$$其中$u$为弦的位移，$x$为弦的位置，$\omega$为本征频率。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
MATLAB语言与应用
12
2019年7月29日10时16分
MATLAB语言与应用
13
【例7-8】
>> y0=lambertw(xx); plot(xx, y0)
2019年7月29日10时16分
MATLAB语言与应用
14
【例7-9】
2019年7月29日10时16分
MATLAB语言与应用
15
2019年7月29日10时16分
MATLAB语言与应用
55
2019年7月29日10时16分
MATLAB语言与应用
56
2019年7月29日10时16分
MATLAB语言与应用
57
本章内容简介
函数名 solve（） fsolve（） optimset（） dsolve（） ode45（） odeset（）
2019年7月29日10时16分
MATLAB语言与应用
42
2019年7月29日10时16分
MATLAB语言与应用
43
【例7-15】
2019年7月29日10时16分
MATLAB语言与应用
44
2019年7月29日10时16分
MATLAB语言与应用
45
附加参数的微分方程求解
【例7-16】
2019年7月29日10时16分
MATLAB语言与应用
22
7.2.1.2 微分方程的解析解方法
2019年7月29日10时16分
MATLAB语言与应用
23
【例7-11】
2019年7月29日10时16分
MATLAB语言与应用
24
2019年7月29日10时16分
MATLAB语言与应用
25
2019年7月29日10时16分
MATLAB语言与应用
MATLAB语言与应用
30
7.2.2 微分方程问题的 数值解法
微分方程问题算法概述 一阶微分方程组的数值解 微分方程转换
2019年7月29日10时16分
MATLAB语言与应用
31
7.2.2.1 微分方程问题算法概 述
2019年7月29日10时16分
MATLAB语言与应用
32
2019年7月29日10时16分
4
7.1.1.2 二元方程的图解法
【例7-2】
2019年7月29日10时16分
MATLAB语言与应用
5
7.1.2 多项式型方程的准解析解法
【例7-3】
2019年7月29日10时16分
MATLAB语言与应用
6
2019年7月29日10时16分
MATLAB语言与应用
7
【例7-4】
2019年7月29日10时16分
2019年7月29日10时16分
MATLAB语言与应用
3
7.1.1 代数方程的图解法
7.1.1.1 一元方程的图解法
【例7-1】
t=3.5203; vpa(exp(-3*t)*sin(4*t+2)+…
4*exp(-0.5*t)*cos(2*t)-0.5)
2019年7月29日10时16分
MATLAB语言与应用
MATLAB语言与应用
46
2019年7月29日10时16分
MATLAB语言与应用
47
2019年7月29日10时16分
MATLAB语言与应用
48
2019年7月29日10时16分
MATLAB语言与应用
49
7.2.2.3 微分方程转换
单个高阶常微分方程处理方法
2019年7月29日10时16分
MATLAB语言与应用
MATLAB语言与应用
8
【例7-5】
2019年7月29日10时16分
MATLAB语言与应用
9
【例7-6】
2019年7月29日10时16分
MATLAB语言与应用
10
【例7-7】
2019年7月29日10时16分
MATLAB语言与应用
11
7.1.3 一般非线性方程数值解
2019年7月29日10时16分
MATLAB语言与应用
16
【例7-10】
2019年7月29日10时16分
MATLAB语言与应用
17
2019年7月29日10时16分
MATLAB语言与应用
18
7.2 微分方程的求解
常系数线性微分方程的解析解方法 微分方程问题的数值解法
2019年7月29日10时16分
MATLAB语言与应用
19
7.2.1 常系数线性微分方程 的解析解方法
线性常系数微分方程解析解的数学描述
微分方程的解析解方法
2019年7月29日10时16分
MATLAB语言与应用
20
7.2.1.1 线性常系数微分方程 解析解的数学描述
2019年7月29日10时16分
MATLAB语言与应用
21
2019年7月29Байду номын сангаас10时16分
50
2019年7月29日10时16分
MATLAB语言与应用
51
【例7-17】
2019年7月29日10时16分
MATLAB语言与应用
52
高阶常微分方程组的变换方法
2019年7月29日10时16分
MATLAB语言与应用
53
【例7-18】
2019年7月29日10时16分
MATLAB语言与应用
54
2019年7月29日10时16分
MATLAB语言与应用
33
2019年7月29日10时16分
MATLAB语言与应用
34
2019年7月29日10时16分
MATLAB语言与应用
35
2019年7月29日10时16分
MATLAB语言与应用
36
2019年7月29日10时16分
MATLAB语言与应用
37
7.2.2.2 一阶微分方程组的数值解
26
处理系数后方程的解：
2019年7月29日10时16分
MATLAB语言与应用
27
2019年7月29日10时16分
MATLAB语言与应用
28
【例7-12】
2019年7月29日10时16分
MATLAB语言与应用
29
7.2.1.3 特殊非线性微分方程的解析解
【例7-13】
2019年7月29日10时16分
四阶五级Runge-Kutta-Felhberg算法
2019年7月29日10时16分
MATLAB语言与应用
38
求解函数
2019年7月29日10时16分
MATLAB语言与应用
39
2019年7月29日10时16分
MATLAB语言与应用
40
2019年7月29日10时16分
MATLAB语言与应用
41
【例7-14】
第7章 代数方程和微分方程的 计算机求解
现代设计与分析研究所
王雷
2019年7月29日10时16分
MATLAB语言与应用
1
主要内容
代数方程的求解 微分方程的求解
2019年7月29日10时16分
MATLAB语言与应用
2
7.1 代数方程的求解
代数方程的图解法 多项式型方程的准解析解法 一般非线性方程数值解