校内培训微分方程求解及MATLAB编程PPT课件.pptx

微分方程求解及MATLAB编程 王义龙 提纲 微分方程模型常微分方程模型及其求解偏微分方程模型及其求解 在许多实际问题中 当直接导出变量之间的函数关系较为困难 但导出包含未知函数的导数或微分的关系式较为容易时 可用建立微分方程模型的方法来研究该问题 在自然界以及工程技术领域中 微分方程模型是大量存在的 它甚至可以渗透到人口问题以及商业预测等领域中去 其影响是广泛的 随时间 空间 变化的数量关系微分方程 含有未知函数的导数 或微分 的方程例 人口模型 种群竞争模型 一 微分方程模型 4 微分方程模型的作用 描述对象特征随时间 空间 的演变过程 分析对象特征的变化规律 预报对象特征的未来性态 研究控制对象特征的手段 根据函数及其变化率之间的关系确定函数 微分方程建模的方法 根据建模目的和问题分析作出简化假设 按照内在规律或用类比法建立微分方程 一 微分方程模型 微分方程建模的对象涉及 改变 变化 增加 减少 衰变 边际 速度 运动 追赶 逃跑 等等词语的确定性连续问题 微分方程建模的基本手段微元法物理定律等 一 微分方程模型 微分方程建模类型 1 常微分方程未知函数为一元函数如 2 偏微分方程未知函数为多元函数如 一 微分方程模型 微分方程模型的一般步骤 一 微分方程模型 一 模型建立1 确定未知函数2 基于机理分析 微元法等 建立关于未知函数的微分方程3 建立定解条件初始条件或边界条件4 写出模型 方程 条件 二 模型求解1 将方程离散化差分格式或有限元 2 将定解条件离散化3 建立问题 模型 的离散化代数方程 差分方程 通常是非线性方程组4 参数的确定5 求解非线性代数方程的算法6 结果7 分析 截断误差 计算误差 灵敏度分析 稳定性 收敛性等 一 微分方程模型 从上面的一般步骤可以看出 微分方程模型的几个要点 1 必须有定解条件 没有定解条件无法求出解 2 方程中参数的确定 通常用最小二乘法 3 对微分方程的求解有两个途径 简单点的可以用Matlab中的ODEPDEPEPEDTOOL等求解 复杂和特殊的方程或区域就需要自己转化为代数方程组求解 此时 非线性方程组的求解算法就是必须的 一 微分方程模型 对参数的估计一般采取两种方法 在没有测试数据的时候查阅文献直接取值 在有测试数据的时候用最小二乘法 一 微分方程模型 比如离散后的非线性方程组为 给出的测试数据是 利用最小二乘可以求出参数 微分方程离散化 差分格式 一阶向前差商 一阶向后差商 一阶中心差商 一 微分方程模型 二阶中心差商 比如热传导方程 用一阶向前 向后 中心差商 二阶中心差商可以得到差分格式为 截断误差为 一 微分方程模型 一 微分方程模型 截断误差为 截断误差为 二 常微分方程模型及其求解 一般的 常 微分方程或微分方程组可以写成 初值问题 定量 定性 解有初等函数表达式 数值解 计算机求解 稳定性 二 常微分方程模型及其求解 1 问题的提出 人口问题是当今世界上最令人关注的问题之一 一些发展中国家的人口出生率过高 越来越威胁着人类的正常生活 有些发达国家的自然增长率趋于零 甚至变为负数 造成劳动力紧缺 也是不容忽视的问题 另外 在科学技术和生产力飞速发展的推动下 世界人口以空前的规模增长 统计数据显示 二 常微分方程模型及其求解 16 我国是世界第一人口大国 地球上每九个人中就有二个中国人 在20世纪的一段时间内我国人口的增长速度过快 如下表 认识人口数量的变化规律 建立人口模型 作出较准确的预报 是有效控制人口增长的前提 下面介绍两个最基本的人口模型 Malthus模型 Logistic模型 二 常微分方程模型及其求解 17 18 2 模型1 Malthus模型 18世纪末 英国人Malthus在研究了百余年的人口统计资料后认为 在人口自然增长的过程中 净相对增长率 出生率减去死亡率为净增长率 是常数 二 常微分方程模型及其求解 二 常微分方程模型及其求解 二 常微分方程模型及其求解 二 常微分方程模型及其求解 二 常微分方程模型及其求解 Matlab求解 利用函数dsolve 命令如下symsN t rN0dsolve DN r N N 0 N0 ans N0 exp r t 二 常微分方程模型及其求解 r 0 2743 10年 xm 4 188 数据拟合 二 常微分方程模型及其求解 曲线拟合可直接利用曲线拟合cftool工具箱 curvefitting 指数增长模型的应用及局限性 与19世纪以前欧洲一些地区人口统计数据吻合 适用于19世纪后迁往加拿大的欧洲移民后代 可用于短期人口增长预测 不符合19世纪后多数地区人口增长规律 不能预测较长期的人口增长过程 19世纪后人口数据 二 常微分方程模型及其求解 分析表明 以上这些现象的主要原因是随着人口的增长 自然资源 环境条件等因素对人口增长的限制作用越来越显著 人口较少时 人口的自然增长率基本上是常数 而当人口增加到一定数量以后 这个增长率就要随着人口的增加而减少 因此 我们将对指数模型关于净相对增长率是常数的基本假设进行修改 2 5模型修改 29 2 利用Matlab直接计算得 symsN t N0Nmr dsolve DN r 1 N Nm N N 0 N0 ans Nm exp Nm log N0 Nm N0 Nm r t Nm 1 化简得 34 x t S形曲线 x增加先快后慢 这个模型称为阻滞增长模型 Logistic模型 Logistic曲线 35 参数估计 用指数增长模型或阻滞增长模型作人口预报 必须先估计模型参数r或r xm 利用统计数据用最小二乘法作拟合 例 美国人口数据 单位 百万 阻滞增长模型 Logistic模型 模型检验 用模型计算2000年美国人口 与实际数据比较 实际为281 4 百万 模型应用 预报美国2010年的人口 加入2000年人口数据后重新估计模型参数 阻滞增长模型 Logistic模型 x 2010 306 0 3 06188亿 2009年 世界国家和地区第3名 次于中国 印度 人口密度31人 平方公里 世界国家和地区第177名 4 传染病模型 问题 描述传染病的传播过程 分析受感染人数的变化规律 预报传染病高潮到来的时刻 预防传染病蔓延的手段 按照传播过程的一般规律 用机理分析方法建立模型 已感染人数 infective i t 每个病人每天有效接触 足以使人致病 人数为 模型1 假设 若有效接触的是病人 则不能使病人数增加 建模 模型2 区分已感染者 infective 和未感染者 susceptible 假设 1 总人数N不变 病人和健康人的比例分别为 2 每个病人每天有效接触人数为 且使接触的健康人致病 建模 日接触率 SI模型 模型2 tm 传染病高潮到来时刻 日接触率 tm t tm di dt最大 卫生水平 模型3 传染病无免疫性 病人治愈成为健康人 健康人可再次被感染 增加假设 SIS模型 3 病人每天治愈的比例为 日治愈率 建模 日接触率 1 感染期 一个感染期内每个病人的有效接触人数 称为接触数 模型3 接触数 1 阈值 感染期内有效接触感染的健康者人数不超过病人数 模型2 SI模型 如何看作模型3 SIS模型 的特例 45 模型4 传染病有免疫性 病人治愈后即移出感染系统 称移出者 Removed SIR模型 假设 1 总人数N不变 病人 健康人和移出者的比例分别为 2 病人的日接触率 日治愈率 接触数 建模 需建立的两个方程 对移出者 模型4 SIR模型 编SIR模型的代码 保存为sir m 根据ode45 函数对微分方程的形式要求 dy sir t y dy必须是列向量 y 1 dy1 dt y 2 dy2 dt y 3 dy3 dt 故sir m第3行意思就是dy 1 dt lamda y1 y2 miu y1 dy 2 dt lamda y1 y2 dy 3 dt miu y1 正是sir模型方程 y1 i y2 s y3 r lamda和miu的值可根据自己的情况设定 functiondy sir t x lamda 0 75 miu 0 25 dy zeros 3 1 dy 1 lamda y 1 y 2 miu y 1 dy 2 lamda y 1 y 2 dy 3 miu y 1 再编写一个运行的m文件 可命名为sirrun m 代码如下 ode45是求解常微分方程最常的函数 其中 sir是sir函数的句柄 0 40 是t的范围 0 040 950 01 分别是I S和R的初始值 t x ode45 sir 0 40 0 040 950 01 t x plot t y 1 t y 2 t y 3 2020 2 8 52 5 用Matlab软件求常微分方程的数值解 t x solver f ts x0 options 1 在解n个未知函数的方程组时 x0和x均为n维向量 m 函数文件中的待解方程组应以x的分量形式写成 2 使用Matlab软件求数值解时 高阶微分方程必须等价地变换成一阶微分方程组 注意 55 注意1 1 建立M文件函数functionxdot fun t x y xdot f1 t x t y t f2 t x t y t 2 数值计算 执行以下命令 t x y ode23 fun t0 tf x0 y0 注意 执行命令不能写在M函数文件中 56 例如 令 注意2 functionxdot fun1 t x y fun1 m xdot f t x t y t x t t x y ode23 fun1 t0 tf x0 y0 M 文件函数形式 y t 是原方程的解 x t 只是中间变量 如果方程形式为 z f t z z 57 解 令y1 x y2 y1 1 建立m 文件vdp1000 m如下 functiondy vdp1000 t y dy zeros 2 1 dy 1 y 2 dy 2 1000 1 y 1 2 y 2 y 1 2 取t0 0 tf 3000 输入命令 T Y ode15s vdp1000 03000 20 plot T Y 1 3 结果如图 例1 解1 建立m 文件rigid m如下 functiondy 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 取t0 0 tf 12 输入命令 T Y ode45 rigid 012 011 plot T Y 1 T Y 2 T Y 3 3 结果如图 图中 y1的图形为实线 y2的图形为 线 y3的图形为 线 例2解微分方程组 59 该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系 没有考虑食饵和捕食者自身的阻滞作用 是Volterra提出的最简单的模型 例3 捕食系统的Volterra方程 m 文件shier m如下 functiondx shier t x dx zeros 2 1 dx 1 x 1 1 0 1 x 2 dx 2 x 2 0 5 0 02 x 1 主程序shark m如下 t x ode45 shier 015 252 plot t x 1 t x 2 plot x 1 x 2 求解结果 左图反映了x1 t 与x2 t 的关系 可以猜测 x1 t 与x2 t 都是周期函数 模型 二 考虑人工捕获 设表示捕获能力的系数为e 相当于食饵的自然增长率由r1降为r1 e 捕食者的死亡率由r2增为r2 e 设战前捕获能力系数e 0 3 战争中降为e 0 1 则战前与战争中的模型分别为 模型求解 1 分别用m 文件shier m和shier m定义上述两个方程 2 建立主程序shark m 求解两个方程 数学模型建立的一般方法 确定所研究的物理量 建立适当的坐标系 划出研究小单元 根据物理定律和实验资料写出该单元与邻近单元的相互作用 分析这种相互作用在一个短时间内对所研究物理量的影响 表达为数学式 简化整理 得到方程 三 偏微分方程模型及其求解 热传导动方程 热传导方程的导出和定解条件 一 热传导方程的导出 给定一空间内物体 设其上的点在时刻的温度为 模型 问题 研究温度的运动规律 分析 两个物理定律和一个公式 1 热量守恒定律 2 傅里叶 Fourier 热传导定律 物体在时间dt内流过面积dS的热量dQ与物体温度沿曲面dS法线方向的方向导数成正比 即 温度变化吸收的热量 通过边界流入的热量 热源放出的热量 为热传导系数 3 热量公式 任取物体内一个由光滑闭曲面所围成的区域 研究物体在该区域内热量变化规律 热传导方程的推导 热量守恒定律 区域内各点的温度从时刻的温度改变为时刻的温度所吸收 或放出 的热量 应等于从时刻到时刻这段时间内通过曲面流入 或流出 内的热量和热源提供 或吸收 的热量之和 即 内温度变化所需要的热量 通过曲面流入内的热量 热源提供的热量 下面分别计算这些热量 1 内温度变化所需要的能量 那么包含点的体积微元的温度从变为所需要的热量为 设物体 的比热 单位质量的物体温度改变 所需要的热量 为 密度为 整个内温度变化所需要的能量 2 通过曲面进入内的热量 由傅里叶热传导定律 从到这段时间内通过进入内的热量为 由高斯公式 知 3 热源提供的热量 用表示热源强度 即单位时间内从单位体积内放出的热量 则从到这段时间内内热源所提供的热量为 由热量守恒定律得 由及的任意性知 三维无热源热传导方程 三维有热源的热传导方程 均匀且各向同性物体 即都为常数的物体 其中 称为非齐次项 自由项 通常称 为非齐次的热传导方程 而称 为齐次热传导方程 定解条件 初始条件和边界条件 初始条件 边界条件 1 第一边界条件 Dirichlet边界条件 特别地 时 物体表面保持恒温 2 第二边界条件 Neumann边界条件 特别地 时 表示物体绝热 3 第三边界条件 D N混合边界条件 其中 表示沿边界上的单位外法线方向的方向导数 注 注意第三边界条件的推导 研究物体与周围介质在物体表面上的热交换问题把一个温度变化规律为的物体放入空气介质中 已知与物体表面接触处的空气介质温度为 它与物体表面的温度并不相同 这给出了第三边界条件的提法 热传导试验定律或牛顿定律 从物体流到介质中的热量和两者的温差成正比 其中比例常数称为热交换系数 流过物体表面的流量可以从物质内部 傅里叶定律 和外部介质 牛顿定律 两个方面来确定 或 即得 1 热传导方程不仅仅描述热传导现象 也可以刻画分子 气体的扩散等 也称扩散方程 注 2 除了三维热传导方程外 物理上 温度的分布在同一个界面上是相同的 可得一维热传导方程 而对于薄片的热传导 可得二维热传导方程 拉普拉斯方程 当我们研究物理中的各类现象 如振动 热传导 扩散等的稳定过程时 由于表达该物理过程的物理量不随时间变化而变化 因此 如果我们考虑的是一个稳定的热场 则可以得到不随时间变化而变化的温度所满足的方程 方程 称为三维拉普拉斯 Laplace 方程或者调和方程 它通常表示成为或者的形式 拉普拉斯方程和泊松方程不仅描述稳定状态下温度的分布规律 而且也描述稳定的浓度分布及静电场的电位分布等物理现象 下面介绍利用有限差分方法来数值求解椭圆型偏微分方程 抛物型偏微分方程 Helmholtz方程是椭圆型偏微分方程中的一种 它满足 椭圆型偏微分方程数值解 自变量取值范围D为 边界条件为 如果g x y 0 则称上述方程为Possion方程 如果同时满足g x y 0和f x y 0 则成为Laplace方程 1 求解原理描述对区域D进行分割 将 范围内的x轴等分 为Mx段 每段长为 同理将y轴分割成My段 每段长为 利用中心差分估计二阶导数 其中 对于自变量区间上的每一个内点 可得到如下等式 其中 在区域D上共有 My 1 Mx 1 个内点 因此可以建立 My 1 Mx 1 个方程 但当Mx My较大时 直接求解 不现实 迭代法可以解决这个问题 迭代法首先把边界条件变形为如下形式 其中 function u x y Helmholtz f g bx0 bxf by0 byf D Mx My MinErr MaxIter 解方程 u xx u yy g x y u f x y 自变量取值区域 D x0 xf y0 yf x y x0 x xf y0 y yf 边界条件 u x0 y bx0 y u xf y bxf y u x y0 by0 x u x yf byf x x轴均分为Mx段 y轴均分为My段 tol误差因子 MaxIter 最大迭代次数 2 求解程序Matlab的实现 x0 D 1 xf D 2 y0 D 3 yf D 4 dx xf x0 Mx x x0 0 Mx dx 构造内点数组dy yf y0 My y y0 0 My dy Mx1 Mx 1 My1 My 1 边界条件form 1 My1u m 1Mx1 bx0 y m bxf y m 左右边界endforn 1 Mx1u 1My1 n by0 x n byf x n 上下边界end 边界平均值作迭代初值sum of bv sum sum u 2 My 1Mx1 u 1My1 2 Mx u 2 My 2 Mx sum of bv 2 Mx My 2 fori 1 Myforj 1 MxF i j f x j y i G i j g x j y i endend dx2 dx dx dy2 dy dy dxy2 2 dx2 dy2 rx dx2 dxy2 ry dy2 dxy2 rxy rx dy2 foritr 1 MaxIterforj 2 Mxfori 2 Myu i j ry u i j 1 u i j 1 rx u i 1 j u i 1 j rxy G i j u i j F i j 迭代公式endendifitr 1end 3 实例分析 functionex1401 f inline x 2 y 2 x y g inline sqrt x x y x0 0 xf 4 y0 0 yf 4 自变量取值范围Mx 30 My 30 等分段数dbx inline x 2 y 2 x y bx0 inline y 2 y 边界条件bxf inline 4 2 cos y y by0 inline x 2 x byf inline 4 2 cos x x D x0 xfy0yf MaxIter 100 MinErr 1e 4 U x y Helmholtz f g bx0 bxf by0 byf D Mx My MinErr MaxIter clf mesh x y U xlabel x ylabel y zlabel U 若g x y 0 则将程序ex1401 m中的代码 g inline sqrt x x y 改为 g inline 0 x y 若g x y 0 f x y 0 则将程序ex1401 m中的代码 g inline sqrt x x y 改为 g inline 0 x y f inline x 2 y 2 x y 改为 f inline 0 x y 抛物型偏微分