范德波方程简介ppt课件.ppt

VanDerPol方程及应用简介 国际化软件人才实验班 Contents 电子科技大学软件学院 总结 Vanderpol方程概述 Vanderpol方程的推导 Vanderpol方程的数值解及分析 Vanderpol方程在生物医学工程中的应用 Vanderpol方程概述 电子科技大学软件学院 BalthazarvanderPol 1889 1959 电子科技大学软件学院 Vanderpol方程概述 电子科技大学软件学院 Vanderpol方程概述 自上世纪20年代vanderpol方程问世以来 它便成了自激极限环振荡系统的原形 已有不少学者对受迫的vanderpol振荡系统的行为进行研究 其中 能够引起科学家较大兴趣的方面是 在外界周期力的作用下vanderpol振荡系统所产生的现象 特别是同步现象以及动力学中的混沌现象 在物理学 生物学 神经学甚至经济学中 vanderpol振荡已经成为描述振荡过程的一种基础模型 BalthazarvanderPol本人对心脏跳动建立了许多电路模 以此型来对心脏跳动的稳定性进行研究 此后 许多科学家建立了许多不同的模型来研究物理学以及生物学上的一些现象 比如 在生物学中 vanderpol方程已经被运用到口胃神经节的胃磨回路的建模当中 在地震学中 对于地理断层中不同地层之间的相互作用的建模当中 vanderpol方程也起到了很大的作用 具体应用 请参考本ppt的后续章节 电子科技大学软件学院 Vanderpol方程的推导 在此我们利用简化的电感耦合晶体管振荡电路来建立非线性振荡器的简化模型 左图是一振荡电路的交流通路 其中的电流和电压都是交流分量 u是集电极电压 ub是基极电压 i是集电极电流 ib是基极电流 图中还标出了各个电流的正方向 考虑到集电极的变化电压 得到以下两个公式 图 电感耦合晶体管振荡器的交流通路 电子科技大学软件学院 Vanderpol方程的推导 图 电感耦合晶体管振荡器的交流通路 根据基尔霍夫定律可得 另外我们可以得到等式 基极电压由下式表示 电子科技大学软件学院 Vanderpol方程的推导 图 电感耦合晶体管振荡器的交流通路 其中a b均为常数 电子科技大学软件学院 Vanderpol方程的推导 图 电感耦合晶体管振荡器的交流通路 其中 这就是著名的vanderpol方程 将某些变量取特殊值后 就变换成了我们最常见的形式 电子科技大学软件学院 Vanderpol的数值解及分析 在本节 将对vanderpol的数值解法做具体的介绍 并对不同的方法的优劣进行适当的分析 以下介绍以下求解常微分方程的常用的单步法中的三种 显式Euler法改进的Euler方法四阶Runge Kutta法 在介绍这三种方法前 我们先利用matlab的simulink对vanderpol方程进行一次仿真 先来领略一下vanderpol方程的函数图像 Vanderpol方程的数值解及分析 matlab仿真 电子科技大学软件学院 图 vanderpol方程的仿真结构图 Vanderpol方程的数值解及分析 matlab仿真 电子科技大学软件学院 图 vanderpol方程的相平面图 图 vanderpol方程y t函数图 Vanderpol方程的数值解及分析 解vanderpol的基本原理 电子科技大学软件学院 1 列出vanderpol方程 降为一阶 初始条件 Vanderpol方程的数值解及分析 解vanderpol的基本原理 电子科技大学软件学院 3 用matlab中的inline 函数来表示vanderpol一阶微分方程组 以下是函数运行结果 ydot inline y 2 y 1 2 1 y 2 y 1 t y ydotydot Inlinefunction ydot t y y 2 y 1 2 1 y 2 y 1 argnames ydot ans t y Vanderpol方程的数值解及分析 显式Euler法 电子科技大学软件学院 每当提到ODE 常微分方程 的数值解 总免不了经典的 欧拉法 在此我们也不例外的先讨论此方法 在此 我们不加推导的给出Euler方程 设 1 Vanderpol方程的数值解及分析 显式Euler法 电子科技大学软件学院 由欧拉公式可得 实际上 y xn yn yn也有误差 它对yn 1的误差也有影响 但这里不考虑此误差的影响 仅考虑方法或公式本身带来的误差 由 2 1 两式相减得欧拉公式的局部截断误差为 故欧拉方法是一阶方法或具有一阶精度 Vanderpol方程的数值解及分析 显式Euler法 电子科技大学软件学院 图 matlabEuler方法函数的活动图 function t y ode Euler f tspan y0 N ifnargin 4 N 0 N100 endifnargin 3 y0 0 endh tspan 2 tspan 1 N t tspan 1 0 N h y 1 y0 fork 1 Ny k 1 y k h feval f t k y k end 使用matlab脚本语言编写Euler法函数 f 要求解的一阶微分方程 组 tspan 函数自变量区间y0 函数初值N 迭代次数 Vanderpol方程的数值解及分析 显式Euler法 电子科技大学软件学院 function t y ode Euler f tspan y0 N ifnargin 4 N 0 N100 endifnargin 3 y0 0 endh tspan 2 tspan 1 N t tspan 1 0 N h y 1 y0 fork 1 Ny k 1 y k h feval f t k y k end 代码的大体框架解释 对输入的参量进行检查 计算步长 并得出自变量向量 赋初值 因为是方程组 采用向量方式赋值 注 由于代码结构非常简单 改进的Euler法和四阶Runge Kutta法和此结构相似 所以后两种方法的代码结构不再赘述 Vanderpol方程的数值解及分析 显式欧拉法 电子科技大学软件学院 调用ode euler函数 ydot inline y 2 y 1 2 1 y 2 y 1 t y tspan 0 20 区间设置y0 0 25 0 5 N 2000 t y ode euler ydot tspan y0 N subplot 121 plot t y 1 xlabel t ylabel y grid 函数y的图像subplot 122 plot y 2 y 1 xlabel y ylabel y grid vanderpol的相平面图tic t y ode euler ydot tspan y0 N time euler toc 计算欧拉方法的运行时间disp time euler Vanderpol方程的数值解及分析 显式欧拉法 电子科技大学软件学院 运行结果分析 测试次数20次 得到Euler方法的运行时间 单位 s 0 82230 73710 74910 74310 75810 76200 73060 74470 72200 76090 75860 73930 76240 77410 77320 71740 75250 76990 79380 7440 平均执行时间为0 7558 Vanderpol方程的数值解及分析 显式欧拉法 电子科技大学软件学院 图 使用Euler方法绘制出的vanderpol方程图像 Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 在介绍改进的Euler法之前 先简单的介绍一下梯形法 下面给出梯形法的公式 上式成为 改进的Euler方法 Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 改进Euler法的局部误差估计 由改进Euler法可得 Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 以上两式相减可以得到 改进的Euler方法的局部截断误差为 故改进的Euler方法是二阶方法 比Euler方法的精度要高 同时 此方法是显示方法 避免了梯形法的可能带来的过量计算 效率较低的缺陷 改进的Euler方法又称欧拉预估 校验方法 具体地说就是先用欧拉公式求出一个近似值 然后利用梯形公式对它进行校正 按照上述思想 可以得到其他预估 校正方法 预估时可以采用中矩公式 抛物线公式进行预估 在此不再赘述 Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 使用matlab编写改进的Euler方法 function t y ode improvedEuler f tspan y0 N ifnargin 4 N 0 N100 endifnargin 3 y0 0 endh tspan 2 tspan 1 N t tspan 1 0 N h y 1 y0 fork 1 Nfk feval f t k y k y k 1 y k h fk y k 1 y k h 2 fk feval f t k 1 y k 1 end f 要求解的一阶微分方程 组 tspan 函数自变量区间y0 函数初值N 迭代次数 图 matlab改进Euler方法函数活动图 Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 调用ode improvedeuler函数 ydot inline y 2 y 1 2 1 y 2 y 1 t y tspan 0 20 y0 0 25 0 5 N 2000 t y ode improvedEuler ydot tspan y0 N subplot 121 plot t y 1 xlabel t ylabel y gridsubplot 122 plot y 2 y 1 xlabel y ylabel y gridtic t y ode euler ydot tspan y0 N time euler toc disp time euler Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 运行结果分析 测试次数20次 得到改进Euler方法的运行时间 单位 s 0 74680 77610 77510 74150 75340 77440 74220 74480 74280 72160 75740 71890 75680 74770 75060 74650 73730 75650 74530 7308 平均执行时间为0 7483 Vanderpol方程的数值解及分析 改进的Euler法 电子科技大学软件学院 图 使用改进Euler方法绘制出的vanderpol方程图像 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 Runge Kuttaf方法的基本思想 Runge Kutta相比之前介绍的两种单步方法来说具有较高的精度 此方法不是通过求导数的方法构造近似公式 而是通过计算不同点的函数值 并对这些函数值作线性组合 构造近似公式 再把近似公式与解的Taylor展开式进行比较 使前面的若干项相同 从而使近似公式达到一定的阶数 在此回顾一下前面介绍的两种方法 Euler法和改进的Euler方法 对于Euler法 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 对于改进的Euler方法 中 令 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 其中 以上就是Runge Kutta方法的基本思想 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 在实际应用中 最常用的是四阶R K法 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 使用matlab编写改进的四阶R K方法 f 要求解的一阶微分方程 组 tspan 函数自变量区间y0 函数初值N 迭代次数 function t y ode RK4 f tspan y0 N varargin ifnargin 4 N 0 N 100 endifnargin 3 y0 0 endy 1 y0 h tspan 2 tspan 1 N t tspan 1 0 N h fork 1 Nf1 h feval f t k y k varargin f1 f1 f2 h feval f t k h 2 y k f1 2 varargin f2 f2 f3 h feval f t k h 2 y k f2 2 varargin f3 f3 f4 h feval f t k h y k f3 varargin f4 f4 y k 1 y k f1 2 f2 f3 f4 6 end Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 图 matlab四阶R K方法函数活动图 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 调用ode RK4函数 ydot inline y 2 y 1 2 1 y 2 y 1 t y tspan 0 20 y0 0 25 0 5 N 2000 t y ode RK4 ydot tspan y0 N subplot 121 plot t y 1 xlabel t ylabel y gridsubplot 122 plot y 2 y 1 xlabel y ylabel y gridtic t y ode euler ydot tspan y0 N time euler toc disp time euler Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 运行结果分析 测试次数20次 得到改进Euler方法的运行时间 单位 s 0 83060 80730 84880 86570 84900 76390 77300 81080 76530 77000 74490 76040 75610 76730 77090 78720 76010 89110 79820 7593 平均执行时间为0 7940 由于Euler 改进的Euler以及RK4时间的复杂度均为O N 所以运行时间差异并不大 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 图 使用改进四阶R K方法绘制出的vanderpol方程图像 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 使用c 实现四阶R K方法 程序的基本思想和流程依照matlab实现的ode RK4方法 将vanderpol方程转换成一次常微分方程组 float f float y floatydot 2 ydot 0 y 1 ydot 1 1 y 1 y 0 y 0 1 y 0 ydot 2 y 2 y 1 2 1 y 1 returnydot Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 a 左边界b 右边界y0 初始值N 迭代次数 voidRunge Kutta floata floatb float y0 intN floatx a x表示当前位置 floaty 2 float K1 NULL K2 NULL K3 NULL K4 NULL temp NULL floath b a N h表示步长 y 0 y0 0 y 1 y0 1 Runge Kutta法的实现函数如下 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 inti j temp float malloc sizeof float 2 K1 float malloc sizeof float 2 K2 float malloc sizeof float 2 K3 float malloc sizeof float 2 K4 float malloc sizeof float 2 cout x 0 x t y y 0 y 1 endl for i 1 i N i for j 0 j 2 j K1 j f y j temp j y j 0 5 h K1 j for j 0 j 2 j K2 j f temp j temp j y j 0 5 h K2 j Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 for j 0 j 2 j K3 j f temp j temp j y j 0 5 h K3 j for j 0 j 2 j K4 j f temp j y j y j h K1 j 2 K2 j 2 K3 j K4 j 6 x a i h cout x i x y y 0 y 1 endl Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 intmain floaty0 2 0 25 0 方程组的初值Runge Kutta 0 20 y0 1000 主函数函数如下 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 程序部分运行结果如下 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 以下我们将对比c 程序的数值结果与matlab计算的数值结果 在此我们主要来观察两个函数图像的拟合程度 为了使两种方法的数值结果具有可比性 在此我们统一以下参数 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 左图蓝色线表示利用c 程序的数值结果绘制的函数图像 红线为matlab的计算结果 误差原因分析 C 程序中 数值结果采用的是float类型 有效位数只有6位 matlab中默认类型为double类型 有效位数可达15位 C 结果精度不够 导致误差 可以看到两结果的曲线拟合的还是很好的 但是 随着t的增大 由于误差的累积两曲线出现了一定的误差 Vanderpol方程的数值解及分析 四阶Runge Kutta法 电子科技大学软件学院 Vanderpol方程在生物医学工程中的应用 电子科技大学软件学院 在本节 将介绍两个关于vanderpol方程的应用 耦合vanderpol振子在心跳建模中的应用Bonhoeffer vanderpol模型在神经细胞薄膜建模中的应用 Vanderpol方程在生物医学工程中的应用 心跳建模 电子科技大学软件学院 将心脏视为耦合的非线性振荡系统这一概念最早追溯到1928年 当时VanderPol和VanderMark使用了三个耦合的松弛振荡 relaxationoscillations 电子电路系统来模拟心脏跳动 VanderPol方程首先描述了在电子电路中的松弛振荡问题 此后 VanderPol方程被经常运用在心脏功能建模的研究中 在心跳模型中 VanderPol方程是一个唯像模型 phenomenologicalmodel 它很好的展示了许多生物学上的特征 例如复杂的周期性 混沌行为等 另外 它还具有以下优点 参数简单 这样可以使得我们研究参数对振动周期的作用更加简单 每一个单独的振荡系统 在没有与其他振荡系统耦合时 有一个稳定的极限环 对于不同的参数数值 振荡系统具有不同的频率 以便我们对比频率变化 作为一种唯像模型 VanderPol方程已被证明具有心脏动力学的许多相似性特征 它可以在描述许多动力学上的特征现象 例如锁模现象 分叉以及混沌现象 Vanderpol方程在生物医学工程中的应用 心跳建模 电子科技大学软件学院 正常的心脏节律是由右心室的一个特殊的叫做SA sino atrial 的节点产生的 SA被视为一个标准的起搏器 normalpacemaker 另外 还有一个起搏器叫做AV atrio ventricular 节 Katholi和他的同事们率先利用VanderPol方程解释了由SA和AV产生的节律之间的相互作用关系 同时考虑了SA以及AV之间耦合的反馈循环 最后对这样的周期性系统做出了分析 由SA节产生的脉冲先后经过心室肌肉组织 AV节 然后传导到心房 其中 心室和心房产生兴奋时有一个毫秒级别的时间延迟 称作AV节阻碍作用 AVBlocks 这种不对称性反应到耦合的振荡系统中 可以有两种形式 SA到AV的单向耦合更精确的一种是 SA与AV的双向耦合 但是这种耦合是非对称的 Vanderpol方程在生物医学工程中的应用 心跳建模 电子科技大学软件学院 SA以及AV节的耦合可以用两个耦合的VanderPol方程来表示 以下给出通式 Vanderpol方程在生物医学工程中的应用 心跳建模 电子科技大学软件学院 下面介绍这个模型的实际应用的意义 心跳节律基本是由两个生理起搏器来产生的 SA AV 这两个起搏器必须同步工作才能在较长的一段时间内维持心跳的稳定 研究发现 心律失常和两个起搏 器工作不同步具有密切的联系 同时这种不同步的状态还可能引起其他心脏方面的疾病 人工起搏器通过对SA节的耦合 调节心律的异常 在实验和临床医学上具有重要意义 Vanderpol方程在生物医学工程中的应用 神经细胞薄膜建模 电子科技大学软件学院 目前 已有大量的针对神经轴突的电子生物特性的研究 其中最著名的也许就是Hodgkin和Huxley所做的工作了 他们对神经轴突的兴奋以及传播进行了量化的研究 并且得出了Hodgkin Huxley方程 HH方程 其中有四个变量V m h n HH模型是利用相平面法来描述神经细胞的兴奋与传导 同样的 一种从VanderPol方程演化而来的模型Bonhoeffer vanderpol BVP 模型也能对以上现象最初很好的阐述 BVP模型只有两个状态 而HH模型有四个状态 所以BVP模型反应出的神经纤维的性质可以在相平面图上表示出来 从BVP模型中 我们可以得到一个生理状态图 physiologicalstatediagram 他不仅可以反映神经冲动的波列 impulsetrains 还可以反映HH模型中很多性质 将相平面图分成不同的区域可以来表示不同的神经纤维的状态 比如静息 兴奋等 所以 在很多情况下BVP模型可以帮助我们理解和解释HH模型 Vanderpol方程在生物医学工程中的应用 神经细胞薄膜建模 电子科技大学软件学院 BVP方程如下 BVP模型旨在表现这一类振动的本质性质 而非去拟合针对某一种情况下的神经纤维反映出的物理数量指标 所以 它简化的代数形式 可以让我们更加集中精力去对此类振动的性质的研究 总结 电子科技大学软件学院 未解决的问题本PPT只对一种形式简单的VanderPol方程进行了介绍 对于其他更加复杂的形式没有解释 在VanderPol方程应用一节中 由于本人能力 时间有限 只做了粗略的介绍 其中VanderPol方程如何在这些问题中具体起何种作用 反映出怎样的性质并未介绍 面临的挑战问题主要集中在VanderPol方程的应用这一块 因为本人未能找到综述性质的并且浅显易懂的文章 而是一些对某些问题较为深入探讨的论文 所以 由于牵扯出的大面积的非振动理论的知识以及对某些模型 比如HH等 没有了解 导致文章理解不深刻 另外一个挑战就是数学知识 虽然大一数学取得了不错的分数 但是长时间未用 知识显得有些生疏 对于公式的理解和推导造成了一定困难 总结 电子科技大学软件学院 以上展示的是本人在这段时间内对VanderPol方程的理解与认识 非常感谢在我从对电路分析 模拟电路 数值分析 非线性系统理论 生物医学工程等知识一无所知 到有一个初步整体的认识过程中 给过我指导和帮助的