两步Runge_Kutta方法求解延迟微分方程的GPL_稳定性_英文_
龙格-库塔方法（Runge-Kutta）
龙格-库塔⽅法（Runge-Kutta）龙格-库塔⽅法（Runge-Kutta）3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的⼀般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分⽤不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值⽅法，若⽤显式单步法(3.2.2)当，即数值求积⽤左矩形公式，它就是Euler法(3.1.2)，⽅法只有⼀阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的⼀种近似，计算时要⽤⼆个右端函数f的值，但⽅法是⼆阶精度的.若要得到更⾼阶的公式，则求积分时必须⽤更多的f值，根据数值积分公式，可将(3.2.1)右端积分表⽰为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)⼀样，⽤前⾯已算得的f值表⽰为(3.2.3)，⼀般情况可将(3.2.2)的表⽰为(3.2.4)其中这⾥均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K⽅法.它每步计算r个f值(即)，⽽k由前⾯(i-1)个已算出的表⽰，故公式是显式的.例i如当r=2时，公式可表⽰为(3.2.5) 其中.改进Euler 法(3.1.11)就是⼀个⼆级显式R-K ⽅法.参数取不同的值，可得到不同公式.3.2.2 ⼆、三级显式R-K ⽅法对r=2的显式R-K ⽅法(3.2.5)，要求选择参数，使公式的精度阶p 尽量⾼，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6) 令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代⼊(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯⼀.因r=2，由（3.2.5）式可知，于是有解(3.2.8)它表明使(3.2.5)具有⼆阶的⽅法很多，只要都可得到⼆阶精度R-K⽅法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常⽤的⼆级R-K⽅法，注意⼆级R-K⽅法只能达到⼆阶，⽽不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个⽅程,加上(3.2.7)的三个⽅程，共计六个⽅程求4个待定参数，验证得出是⽆解的.当然r=2,p=2的R-K⽅法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，⼀般不再给出.对r=3的情形，要计算三个k值，即其中将按⼆元函数在处按Taylor公式展开，然后代⼊局部截断误差表达式，可得可得三阶⽅法，其系数共有8个，所应满⾜的⽅程为这是8个未知数6个⽅程的⽅程组，解也是不唯⼀的，通常.⼀种常见的三级三阶R-K⽅法是下⾯的三级Kutta⽅法：(3.2.11)附：R-K 的三级Kutta ⽅法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满⾜c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K ⽅法及步长的⾃动选择利⽤⼆元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K ⽅法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，⽽33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
matlab龙格库塔方法求解二元二阶常微分方程组
matlab龙格库塔方法求解二元二阶常微分方程组文章标题：深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用摘要：在科学与工程领域，常常需要求解复杂的微分方程组，而matlab作为一种强大的数学工具，提供了许多求解微分方程组的方法。
本文将深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用，以便读者全面理解该方法并能灵活应用于实际问题中。
正文：一、介绍龙格库塔方法龙格-库塔法（Runge-Kutta methods）是一种数值求解常微分方程的方法，通过将微分方程的解进行离散化，将微分方程转化为差分方程，从而进行数值求解。
龙格库塔方法通过迭代计算，能够得到微分方程的数值解，广泛应用于科学计算和工程技术领域。
二、matlab中的龙格库塔方法在matlab中，龙格库塔方法通过ode45函数实现，该函数能够对一阶或高阶常微分方程进行数值求解。
用户可以通过设定初始条件、微分方程表达式，以及积分区间等参数，快速得到微分方程的数值解。
ode45函数采用自适应步长的方式进行求解，能够有效解决微分方程解的数值稳定性和精确度问题。
三、龙格库塔方法在求解二元二阶常微分方程组中的应用考虑如下形式的二元二阶常微分方程组：$$\begin{cases}y_1' = f_1(t, y_1, y_2) \\y_2' = f_2(t, y_1, y_2)\end{cases}$$其中，$y_1(t)$和$y_2(t)$是未知函数，$f_1(t, y_1, y_2)$和$f_2(t,y_1, y_2)$分别表示其对应的函数表达式。
通过matlab中的ode45函数，可以将该二元二阶常微分方程组转化为一阶常微分方程组的形式，然后利用龙格库塔方法进行数值求解。
设定初始条件$y_1(0) = y1_0, y_2(0) = y2_0$，对应的一阶方程组为：$$\begin{cases}u_1' = u_3 \\u_2' = u_4 \\u_3' = f_1(t, u_1, u_2) \\u_4' = f_2(t, u_1, u_2)\end{cases}$$其中，$u_1(t) = y_1(t), u_2(t) = y_2(t), u_3(t) = y_1'(t), u_4(t) =y_2'(t)$，通过ode45函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。
Runge-Kutta方法关于无穷延迟系统的稳定性分析
Ａ （ ） ， ６ （ ， ， ｂ ， Ｃ ∑ 口， １ ５ 一口 ； 一 ６ ｂ …，） ｉ ｓ 一 ≤≤ ，
Ｊ １ ＝
其 中每个 Ｃ 是 不等 的． 对于 （ ． ） １ ５ 中延迟项 ， 用 Ｚ ｎ ａｏ的 自然连 续伸 展 的技 巧 ， 项式 分别 １ ４ 和（ ． ） 借 ｅｎｒ 多
对 比系统
， 一 上 ＋ （ ‘ ｙ Ｍ ｙ（ － ｒ）＋ ｔ １ ￡ 一 ￡ ０’ ＞
（．） １ ３
ｌ ￡ （）一 ｇ（）， 一 ｍａ ｖ ≤ ｔ ０， ￡ ｘ／ ≤
系统 （ ． ） １ ３ 因为是常延 迟则 没有 以上 的 问题 出现． 以我 们 在 以 下 的讨 论 中需 要 借 用一 些 技 巧将 系 统 所 （． ） １ ２ 转化为 系统 （ ． ） １ ３ 的形式 ， 以解决 计算 中的存 储 问题 ．
维普资讯
第２ ３卷 第 ４期
２０ ０ ７年 ８月
大 学 数 学
ＣＯ ＬＬＥＧ Ｅ Ａ Ｔ Ｈ ＥＭ ＡＴ Ｉ Ｍ ＣＳ
Ｖｏ ． ３， ． １２ № ４
Ａ ｕ ．２ ０ ｇ ０７
， ●● ●
●，、 ●● ●●ｌ
Ｒｕ ｇ — ｔａ方 于无 穷延 迟 系统 的稳 定 性 分 析 ｎ ｅＫｕｔ 法关
∑Ｋ ，ｔ ， （￡ Ｙ ｑ） ＞０
＝
１
（ ．２ １ ）
其 中 Ｌ， ， Ｍ ， ∈ ｄｄ Ｙ ，ｑ ∈（ ，） ｚ １ ２ … ， ） Ｍ Ｋ， Ｋ ￣ ， ∈ ０ １ （一 ， ， 为常数 ．
研 究表明， （一 Ｍ ） ，Ｌ＜ 时， Ｄ・对应矩 谱半径， ・为谱 标， 当ＰＬ ∑ ＜１ （） ０ 其中Ｉ ） ａ （ 阵的 ａ ） 横坐 系 （
求解延迟微分方程的ROSENBROCK方法的渐近稳定性
求解延迟微分方程的ROSENBROCK方法的渐近稳定性曹学年;刘德贵;李寿佛
【期刊名称】《系统仿真学报》
【年(卷),期】2002(14)3
【摘要】数值求解延迟微分方程的Runge-kutta方法和θ-方法已经有了较深入的研究。
本文适当改造求解常微分方程的Rosenbrock方法,构造了一类求解延迟微分方程的Rosenbrock方法,证明了这类方法是GP-稳定的,而且这类方法的GP-稳定性与求解常微分方程的Rosenbrock方法的A-稳定性等价。
数值试验表明这类方法是有效的。
【总页数】3页(P290-292)
【关键词】延迟微分方程;Rosenbrock方法;渐近稳定性;数值计算
【作者】曹学年;刘德贵;李寿佛
【作者单位】湘潭大学数学系;北京计算机应用和仿真技术研究所
【正文语种】中文
【中图分类】O241.81
【相关文献】
1.多比例延迟微分方程Rosenbrock方法的渐近稳定性 [J], 刘建国;甘四清
2.延迟微分方程指数Rosenbrock方法的渐近稳定性 [J], 王世英;邢慧
3.求解比例延迟微分方程的Rosenbrock方法的渐近稳定性 [J], 刘建国
4.比例延迟微分方程组Rosenbrock方法的渐近稳定性 [J], 刘建国;甘四清
5.求解中立型比例延迟微分方程组Rosenbrock方法的渐近稳定性 [J], 刘建国;甘四清
因版权原因，仅展示原文概要，查看原文内容请购买。
龙格库塔方法及其matlab实现
龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。
MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MatLab是英文MATrix LABoratory（矩阵实验室）的缩写。
在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
关键词：龙格-库塔 matlab 微分方程1.前言1.1：知识背景龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。
该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。
Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。
经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。
从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。
MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。
龙格-库塔(Runge-Kutta)方法
——常微分方程及方程组的数值解 常微分方程及方程组的数值解
基本概念 只有一个自变量的微分方程为常微分方程 ODE, 否则称为偏微分方程 PDE。 。 方程中未知函数导数的最高阶数称为方程的阶
[1] [2] C. Runge, "Ueber die numerische Auflsung von Differentialgleichungen" Math. Ann. , 46 (1895) pp. 167–178 W. Kutta, "Beitrag zur naherungsweisen Integration von Differentialgleichungen" Z. Math. und Phys. , 46 (1901) pp. 435–453
2. 递推化 在具有唯一解的条件下， 在具有唯一解的条件下，通过步进法逐步计算出解 在一系列离散点上的值。从而得到原 的数值近似解。 在一系列离散点上的值。从而得到原ODE的数值近似解。 的数值近似解
初值问题解法 讨论一阶ODE（高阶可化为一阶ODEs） （高阶可化为一阶 讨论一阶 ） 的初值问题。初值问题可以一般地写成： 的初值问题。初值问题可以一般地写成：
欧拉（ 欧拉（Euler）方法 ）
Euler方法是求解初值问题的最简单方法， 方法是求解初值问题的最简单方法， 方法是求解初值问题的最简单方法 精度差。然而对理论分析很有用。 精度差。然而对理论分析很有用。Runge-Kutta 法是对Euler法的改进 法是对 法的改进
Euler方法： 方法： ym+1 = ym + hK1 + O(h ) K1 = f ( xm , ym )
K 1 = f ( xm , ym ) K i = f ( xm + hai , ym + h∑ bil K l )
两步Runge-Kutta方法求解中立型延迟积分微分代数方程的稳定性
两步Runge-Kutta方法求解中立型延迟积分微分代数方程的稳定性王倩;丛玉豪;李建军【摘要】研究了两步 Runge-Kutta方法求解中立型延迟积分微分代数方程的稳定性,并证明了在某些条件下,A稳定的两步Runge-Kutta方法求解中立型延迟积分微分代数方程可以保持它的渐进稳定性.【期刊名称】《上海师范大学学报（自然科学版）》【年(卷),期】2010(039)006【总页数】7页(P551-557)【关键词】渐进稳定性;中立型延迟积分微分代数方程;两步;Runge-Kutta方法【作者】王倩;丛玉豪;李建军【作者单位】上海师范大学,数理学院,上海,200234;上海师范大学,数理学院,上海,200234;上海师范大学,数理学院,上海,200234【正文语种】中文【中图分类】O241.81 IntroductionConsider the d-dimensional delay integro differential equation:My′(t)=f(t,y(t),y(t-τ),y′(t-τ),g(t,s,y(s))ds),(1.1)where f:R+×Rd×Rd×Rd×Rd→Rd is continuous,τ(>0)is a constant lag and the matrix M∈Rd×d may be singular.This system can be found in a wide variety of scientific and engineering fields such as biology,physics,ecology and so on[5-6].Particularly,it is observed that the delay-integro differential-algebraic system(when M is singular)plays an important role in modeling many different phenomena of circuit analysis and chemical process simulation,which have a comprehensive list in[2,9].Recently,there has been a growing interest in developing numerical methods for solving the system(1.1)[1,2,9].Some papers have discussed the stability of Runge-Kutta methods for the linear neutral delay differential-algebraic equations: Au′(t)+Bu(t)+Cu′(t-τ)+Du(t-τ)=0,(1.2)where the matrix A is singular[9] other papers have discussed the stability analysis of neutral delay integro differential algebraic system up to now.Recently Zhao[3] considered the stability of θ-methods and BDF methods for linear neutral delay integro differential equation:Au′(t)+Bu(t)+Cu′(t-τ)+Du(t-τ)+Gu(σ)dσ=0,(1.3)where the matrix A may be singular.It was proved that every linear θ methods with θ∈(1/2,1]and A-stability BDF method preserve the delay-independent stability of its exact solutions in[3].In this paper,we focus on stability of two step Runge-Kutta methods (TSRK) for the linear neutral delay integro differential algebraic equations(NDIDAE).It is proved that thetwo step Runge-Kutta methods preserve the delay-independent stability under some conditions if TSRK is A-stable.2 Asymptotic stability results for NDIDAEs2.1 Asymptotic stability of NDIDAEsIn this subsection,we consider the following linear system of neutral delay differential-algebraic equations of the form(2.1)where A,B,C,D,G∈d×d,A is a singular matrix, τ is a given positive delay constant and τ=mh,m∈N,φ(t)denotes a given vector-valued function,u(t)is a vector-valued unknown function to be solved for t≥0.Forsolvability,which is essentially the existence and uniqueness of its exact solutions,an important conclusion is given by the followingLemma[1,2,9].The characteristic polynomial of the NDIDAEs(2.1) is given as follows:π(z)=det[zA+B+Cze-zτ+De-zτ+Gz-1(1-e-zτ)].(2.2)Since z=0 is a removable singularity of Eq.(2.2),then we can regard π(z)as an entire function withπ(0)=det(B+D+τ G).(2.3)the zero solution of Eq.(2.1)is asymptotic stability if and only if all roots of π(z)=0 have negative real parts,Hence,its delay-independent stability is characterized as(see P101 of [4])π(z)=0⟹R(z)≤-r<0 with r∈R and r>0 forany τ≥0.Lemma 2.1[1] Under conditions:•a(γ)>c(γ) for any γ∈Cdwith a(γ)≠0;•b(γ)>d(γ) for any γ∈Cdand a(γ)=c(γ)=0;the system (2.1)is asymptotically stable if and only if(Ⅰ) det(B+D+τG)≠0,(Ⅱ) det(z2A+zB+G)=0,z≠0⟹R(z)<0,(Ⅲ)where ρ(·)denotes the spectral radius of the matrix and a(γ)=〈γ,Aγ〉,b(γ)=〈γ,Bγ〉,c(γ)=〈γ,Cγ〉,d(γ)=〈γ,Dγ〉.2.2 Asymptotic stability of TSRK methods for NDIDAEsIn this subsection,we will concern the asymptotic stability of TSRK methods for NDIDAEs.2.2.1 Brief introduction to TSRK methodsConsider the TSRK methods of the form(2.4)where τj=tn+h,=tn-1+h and 0≤θ≤1.We will use the following table:where = and (+)=1+θ.In [4],the authors investigate stability propertiesof(2.4)with respect to the basic test equationx′(t)=αx(t), t≥0,Re α<0,and obtainxi+1=R(α,θ)xi+S(α,θ)xi-1,(2.5)whereThe characteristic polynomial of (2.5) ish(z)=z2-R(α,θ)z-S(α,θ).(2.6)The stability region of the TSRK methods (2.4) is the set of all points α for which the roots of h(z)are inside or on the unit circle with those on the unit circle being simple. The method is said to be A-stable if its stability region contains the negative half plane. So if h(z)is a Schur polynomial for any α,Re(α)<0,the TSRK methods is A-stable for ODEs.2.2.2 Asymptotic stability of TSRK methods for NDIDAEsFirst we give a definition of asymptotic stability for numerical methods. Definition 2.1 A numerical method for asymptotically stable neutral system(2.1)is also called asymptotically stable if the numerical solution satisfies un=0.Applying the two-step methods (2.4) to (2.1),we havewhere Un,i and unare approximations to u(tn+cih) and u(tn).Un-m,i is an approximation to u(tn+cih-τ) and Vn,i is an approximation to Gu(σ)dσ.For n=0,1,…,define vnby vn=hGUn-k,j.In view of Eqs.(2.7),we can write (2.7) in the formwhere yn=(,)T,and Zn=(,,…,,)T,=(U,V,…,U,V)T,xn=(,,Z)T then we can express Eq.(2.8-2.10)as(2.11)The characteristic polynomial of (2.11) is given by(2.12)where(2.13)Theorem 2.1 Assume that the TSRK method is A-stable and the real part of any eigenvalue of matrix is non-negative.Then the numerical method is asymptotically stable for system(2.1) with conditions of Lemma 2.1 and det(A+η C)≠0 for 0<η≤1.Proof It follows from Lemma 2.1,det(A+ηC)≠0 and the followin g relation(2.14)f(z,ξ)=(z2A+zB+G)+(z2C+zD-G)ξ.From Lemma 2.1 we haveThenso σ(-(K+Nξ)-1(L+Mξ))⊂C-∪0.From [1,9],we know that under the condition of Lemma 2.1,T33 and K+z-mN is nonsingular for z≥1. Hence p(z)can be rewritten as(2.15)where det(T22)=det(T33)≠0.Let Q(z)=-(K+z-mN)-1(L+z-mM),⊗⊗(Q(z))-1(es⊗I2d),⊗⊗(Q(z))-1(es⊗I2d).By direct calculation,we have(2.16)If the difference equation(2.1) is asymptotically stable,then p(z)is a Schur polynomial.Assume that p(z) is not a Schur polynomial,which implies that there exists a root z of characteristic equatious with z≥1.Because T33 and K+z-mN is nonsingular for z≥1,it is easy to see that p(z)=0⟹h(z)=0.Let J denote Jordan form of matrix Q(z).Then assume λi(i=1,…,s)are the eigenvalues ofQ(z).There fore,there exists at least one λi∈σ(Q(z))⊂C-∪0 such thatIf z>1 then Eq.(2.16)contradicts A-stability of the methods.But if z=1 and R(λi)<0,Eq.(2.16) also contradicts A-stability.Therefore λi=0 which yieldsR(θ)=1-θ,S(θ)=θ,h(z)=1-z-1(1-θ)-z-2θ=0⟹z=1.If 1 is a characteristic root of Eqs.(2.7),then there existu,Ui(i=1,…,s),(i=1,…,s)∈Rd which are not all zero,such that un=u,Un,i=Ui,= satisfy Eqs.(2.7a-2.7f),i.e.From the equations we have vn=vn-1,so from (2.17a-2.17b) we have τG=τGUj.From(2.17c-2.17d),that(Is⊗⊗andBecause σ()⊂ C+∪0 and σ(-h(A+C)-1(B+D))⊂C-∪0.Sodet(Is⊗⊗(B+D))≠0.From Cramer rule we have U-=0.Thus(2.17c)and(2.17e)imply that⊗(A+C)-1(B+D))-1u(es⊗I2d)anduId=(R(θ)+S(θ ))u.Because two step Runge-Kutta methods are A-stability,so det(I2d-R(θ)-S(θ))≠0.Hence u=0 and U=0.But this contradicts theassumption,therefore,1 is not the characteristic root of Eqs.(2.8)-(2.10).Then p(z)is a schur polynomial.References:[1] XU Y,ZHAO J J .Stability of Runge-Kutta methods for neutral delay-integro-differential-algebraic system[J].Mathematics and Computers inSimulation,2008,79(1):571-583.[2] PETZOLD L R.Numerical solutions of differential-algebraicequations[J].Advances in Numerical Analysis Ⅳ,1995,32(5):1635-1657. [3] ZHAO J J,XU Y,LIU M Z.Stability analysis of numerical methods for linear neutral Volterra delay-integro-differential system[J].Appl Math Comput,2005,167(2):1062-1079.[4] KUANG J X,CONG Y H.Stability of Numerical Methods for Delay Differential Equtations[M].Beijing:Science Press,2005.[5] HAIRER E,WANNER G.Solving ordinary differential-algebraic equations Ⅱ stiff and differential-algebraic problems[M].NewYork:Springer,1996. [6] BRUNNER H.The numerical solutions of neutral Volterra integro-differential equations with delay arguments[C].Auckland:Proceedings SCADE’93,1993.[7] GAN S Q,ZHENG W M.Stability analysis of Runge-Kutta methods for delay integro differential equtations[J].TsingHua Science and Technology,2004,9(1):185-188.[8] LIU M Z,SPIJKER M N.The Stability of θ methods in the numeri cal solution[J].IMA Numer Anal,1990,10(1):31-34.[9] HU G,MITSUI T.Stability of numerical methods for systems of neutral delay differential equations[J].BIT,1995,35(4):504-515.[10] REICH S.On the local qualitative behavior of differential-algebraic equation[J].Circuits Systems and Signal Processing,1995,14(4):427-443. [11] ZHU W J,PETZOLD L R.Asymptotic stability of linear delay differential-algebraic equtations and numerical methods[J].Appl NumerMath,1997,24(2):247-264.。
随机微分方程Runge-Kutta方法的矩指数稳定及矩渐近稳定性
随机微分方程Runge-Kutta方法的矩指数稳定及矩渐近稳定性张雨馨【摘要】考虑逼近随机微分方程的1.5阶Runge-Kutta法的矩指数稳定性和矩渐近稳定性,对于标量线性检验方程,证明了随机Runge-Kutta法的矩指数稳定性和矩渐近稳定性是一致的,并给出了这两种稳定性的存在条件.%Moment exponentical and moment asymptotic stabilities of Runge-Kutta method with strong order 1.5 for stochastic differential equations was studied. We have proved that for scalar linear test equations, these stabilities of Runge-Kutta method are equivalent and we gave the conditions for these stabilities existing.【期刊名称】《吉林大学学报（理学版）》【年(卷),期】2012(050)001【总页数】2页(P67-68)【关键词】Runge-Kutta方法;矩指数稳定;矩渐近稳定【作者】张雨馨【作者单位】吉林大学数学学院,长春130012;哈尔滨工程大学理学院,哈尔滨150001【正文语种】中文【中图分类】O211.63对于随机微分方程, 常用的数值方法是Euler-Maruyama格式和Milstein近似, 关于它们的稳定性研究目前已有许多结果[1-6]. 但Euler-Maruyama方法的强收敛阶仅为0.5, Milstein方法的强收敛阶为1. Burrage等[7]构造了1.5阶强收敛的随机Runge-Kutta方法, 并在文献[8]中修正了所构造的算法. 本文研究这类算法的p 阶矩指数稳定性和p阶矩渐近稳定性.Stranovich随机微分方程初值问题为:dy(t)=f(y(t))dt+g(y(t))∘dW(t), y(0)=y0∈Rn, (1)其中: f(y(t))为漂移项; g(y(t))为扩散项; W(t)为标准Wiener过程. 做稳定性分析时, 一般只需考虑标量线性检验方程dy(t)=ay(t)dt+by(t)∘dW(t), (2)其中a,b∈R . 选择h为步长, 用yn表示y(tn)的近似值, 将文献[9]提出的1.5阶Runge-Kutta法应用于随机微分方程(1)上, 有(3)其中:是s×s矩阵; α,γ(1),γ(2)是s维向量; I是单位矩阵; ⊗表示Kronecker积, 而F(Y)=(f(Y1)T,…, f(Ys)T)T, G(Y)=(g(Y1)T,…,g(Ys)T)T, e=(1,…,1).将Runge-Kutta法(3)应用到检验方程(2)上, 可得Y=yne+(haA+b(B(1)J1+B(2)(J10/h)))Y ⟹ Y=(I-haA-b(B(1)J1+B(2)(J10/h)))-1yne,因此,yn+1=yn+hαTaY+(γ(1)TJ1+γ(2)T(J10/h))bY ⟹ yn+1=Rn+1yn, (4)其中Rn+1=1+(haαT+b(γ(1)TJ1+γ(2)T(J10/h)))(I-haA-b(B(1)J1+B(2)(J10/h)))-1e.下面给出矩指数稳定和矩渐近稳定的定义[6,9-10].定义1 给定一个步长h>0, 如果对于所有的初值y0∈Rd, 都有则应用在方程(1)上的数值方法称为p阶矩指数稳定的. 当p=2时, 这种稳定性称为均方指数稳定.定义2 如果对于所有的初值y0∈Rd, 都有则应用在方程(1)上的数值方法称为p 阶矩渐近稳定的. 当p=2时, 这种稳定性称为均方渐近稳定.定理1 对于检验方程(2), Runge-Kutta法(3)的p阶矩指数稳定和p阶矩渐近稳定是一致的, 且充要条件均为E()<1.证明: 由式(4)通过迭代可知=…y0. (5)因为{R1,…,Rn}是独立同分布的随机变量序列, 所以对式(5)的两端同时取数学期望可得E()=(E())ny0, 计算此式的极限可得⟹⟹ E()<1,从而由定义1和定义2可知结论成立.衷心感谢吉林大学数学学院李勇教授的悉心指导.参考文献【相关文献】[1] Buckwar E, Kelly C. Towards a Systematic Linear Stability Analysis of Numerical Methods for Systems of Stochastic Differential Equations [J]. SIAM J Numer Anal, 2010, 48(1): 298-321.[2] Buckwar E, Sickenberger T. A Comparative Linear Mean-Square Stability Analysis of Maruyama and Milstein-Type Methods [J]. Math Comput Simu, 2011, 81: 1110-1127. [3] Higham D J, MAO Xue-rong, YUAN Cheng-gui. Almost Sure and Moment Exponential Stability in the Numerical Simulation of Stochastic Differential Equations [J]. SIAM JNumer Anal, 2007, 45(2): 592-609.[4] Saito Y, Mitsui T. Mean-Square Stability of Numerical Schemes for Stochastic Differential Systems [J]. Vietnam J Math, 2002, 30: 551-560.[5] Saito Y, Mitsui T. Stability Analysis of Numerical Schemes for Stochastic Differential Equations [J]. SIAM J Numer Anal, 1996, 33(6): 2254-2267.[6] Higham D J, MAO Xue-rong, Stuart A M. Exponential Mean-Square Stability of Numerical Solutions to Stochastic Differential Equations [J]. LMS J Comput Math, 2003, 6: 297-313.[7] Burrage K, Burrage P M. High Strong Order Explicit Runge-Kutta Methods for Stochastic Ordinary Differential Equations [J]. Appl Numer Math, 1996, 22(1/2/3): 81-101.[8] Burrage K, Burrage P M. Order Conditions of Stochastic Runge-Kutta Methods by B-Series [J]. SIAM J Numer Anal, 2000, 38(5): 1626-1646.[9] Burrage P M. Runge-Kutta Methods for Stochastic Differential Equations [D]: [Ph D Thesis]. Brisbane, Australia: Department of Mathematics, University of Queensland, 1999.[10] MAO Xue-rong. Stochastic Differential Equations and Applications [M]. Chichester, UK: Horwood, 1997.。
function chap1_rungekutta_method
function chap1_rungekutta_method如何使用龙格-库塔方法（Runge-Kutta method）数值求解常微分方程（Ordinary Differential Equations）。
第一步，我们先来复习一下龙格-库塔方法的原理和基本步骤。
龙格-库塔方法是一种常用的数值求解常微分方程的方法，它基于以下的数值积分思想：通过逐步更新的方式，将微分方程的解近似地表示为离散的数值解。
具体而言，龙格-库塔方法主要包括四个基本步骤：选择步长、计算中间点的斜率、计算权重因子、更新解。
首先，我们需要选择一个适当的步长，用来控制数值求解的精度和效率。
步长越小，数值解越接近真实解，但计算量也越大。
接下来，我们以一个简单的一阶常微分方程为例，介绍如何使用龙格-库塔方法求解。
假设有一个方程dy/dx = f(x, y)，我们需要求解在给定的初始条件下，x = x0时y的值。
第二步，我们需要计算中间点的斜率。
在龙格-库塔方法中，我们通过计算四个斜率来逼近真实值。
具体而言，我们在当前点x处计算斜率k1，然后计算x + 0.5*h处的斜率k2，再计算x + 0.5*h处的斜率k3，最后计算x + h处的斜率k4。
其中，h为步长。
斜率的计算方式为k1 = f(x,y)，k2 = f(x + 0.5*h, y + 0.5*h*k1)，k3 = f(x + 0.5*h, y + 0.5*h*k2)，k4 = f(x + h, y + h*k3)。
第三步，我们需要计算权重因子。
龙格-库塔方法通过加权平均斜率来逼近真实值。
具体而言，我们计算y在下一个步长x + h处的更新值y'，其中y' = y + (1/6)*(k1+2*k2+2*k3+k4)*h。
最后，我们更新解。
将y更新为y'，即y = y'。
以上就是使用龙格-库塔方法求解常微分方程的基本步骤。
接下来，我们通过一个具体的例子来演示如何使用该方法。
中立型多延迟积分微分方程Runge—Kutta方法的散逸性
收 稿 日期 ：２ ０ １ ２ — １ ２ — １ ５
基 金 项 目 ：安徽 省 教 育厅 资 助 项 目 （ＫＪ ２ ０ １ ２ Ｚ３ ６ ７） 作 者 简 介 ：王 素 霞 （１ ９ ８ ２ 一 ） ，女 ，河 南周 口人 ，助 教 。硕 士 ，研 究 方 向 为 微 分 方 程数 值 解．
性 ．本 文 研 究 中 立 型多 延 迟 积 分 微 分 方 程 Ｒ ｕ ｎ ｇ ｅ ． Ｋｕ ｔ ｔ ａ 方 法 的 散逸 性 ，希 望进 一 步完 善 和推 广 已有 的
关 于 中 立 型延 迟 系 统 的 散 逸 性 结 果 ．
１ 中立 型 多 延 迟 积 分 微 分 方 程 的 Ｒｕ ｎ ｇ ｅ — Ｋ ｕ ｔ ｔ ａ方 法
给 出 了 Ｒｕ ｎ ｇ ｅ － Ｋｕ ｔ ｔ ａ方 摘 要 ：研 究 了 中 立 型 多 延 迟 积 分 微 分 方 程 Ｒ ｕ ｎ ｇ ｅ — Ｋ ｕ ｔ ｔ ａ方 法 的 散 逸 ， ，
法的数值散 选性 结果．
关键 词 ： 中 立 型 多 延 迟 积 分 微 分 方 程 ；Ｒｕ ｎ ｇ ｅ － Ｋｕ ｔ ｔ ａ方 法 ；散 逸 性
Ｗ ＡＮＧ Ｓ ｕ －ｘ ｉ ａ ， ＸＵ Ｙｉ ｎ ｇ
（ Ｄ ｅ ｐ ａ ｒ ｔ ｍｅ ｎ ｔ ｏ ｆ Ｍａ ｔ ｈ ｅ ｍａ ｔ ｉ ｃ ｓ ａ ｎ ｄ Ｃ ｏ ｍｐ ｕ ｔ ａ ｔ ｉ ｏ ｎ ａ ｌ Ｓ ｃ ｉ ｅ ｎ ｃ ｅ ， Ｈｕ ａ ｉ ｎ ａ ｎ Ｎｏ ｒ ｍａ ｌ Ｕｎ ｉ ｖ ｅ ｒ ｓ ｉ ｔ ｙ ， Ｈｕ ａ ｉ ｎ ａ ｎ ２ ３ ２ ０ ３ ８ ， Ｃ ｈ ｉ ｎ ａ ）
Ａｂｓ ｔ ｒ ａ ｏ ｔ ： Ｔｈ ｉ ｓ ｐ ａ ｐ ｅ ｒ ｄｅ ａ ｌ ｓ ｗｉ ｔ ｈ ｔ ｈ ｅ ｄ ｉ ｓ ｓ ｉ ｐ ａ ｔ ｉ ｖ ， ｉ ｔ ｙ ｏ ｆ ｔ ｈ ｅ Ｒｕ ｎｇ ｅ — － Ｋｕ ｔ ｔ ａ ｍｅ ｔ ｈ ｏ ｄｓ ｆ ｏ ｒ ｎｅ ｕ ｔ ｒ ａ ｌ ｍｕ ｌ ｔ ｉ — － ｄ ｅ ｌ ａ ｙ ｉ ｎ ｔ ｅ ｇ ｒ ０ 一ｄ ｉ ｆ ｆ ｅ ｒ ｅ ｎ ｔ ｉ ａ ｌ ｅ ｑ ｕ ａ ｔ ｉ ｏ ｎ ｓ ａ ｎ ｄ ｇ ｉ ｖ ｅ ｓ ｔ ｈ ｅ ｄｉ ｓ ｓ ｉ ｐ ａ ｔ ｉ ｖ ｉ ｔ ｙ ｒ ｅ ｓ ｕ ｌ ｔ ｏ ｆ ｔ ｈ ｅ Ｒｕ ｎｇ ｅ —Ｋｕ ｔ ｔ ａ ｍｅ ｔ ｈ ｏ ｄ ｓ ． Ｋｅ ｙ ｗｏ ｒ ｄ ｓ ：ｎ ｅ ｕｔ ｒ ａ ｌ ｍｕｌ ｔ ｉ －ｄ ｅ ｌ ａ ｙ ｉ ｎ ｔ ｅ ｇ ｒ ｏ － ｄ ｉ ｆｅ ｒ ｅ ｎ ｔ ｉ ａ ｌ ｅ ｑｕ ａ ｔ ｉ ｏ ｎ ｓ ； Ｒｕ ｎ ｇ ｅ － Ｋｕ ｔ ｔ ａ ｍｅ ｔ ｈ ｏ ｄ ｓ ；ｄ ｉ ｓ ｓ ｉ ｐ ａ ｔ ｉ ｖ ｉ ｔ ｙ
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
上海师范大学学报 ( 自然科学版 ) Jou rnal of Shanghai N or m al U n iversity( N atural S cien ces)
V ol. 39, N o. 4 A ug. , 2 0 1 0
The GPL stab ility of the two step Runge Kutta m ethods for delay differential system s
[ 7, 8]
: ( 2a ) ( 2b ) ( 2c )
∀
s
aijK n, j ),
s ij n- 1 , j
i= 1 , 2 , #, s, ), i= 1 , 2 , #, s, n = 1 , 2 , #, ! 1 .
j= 1
K n- 1, i = hf ( tn- 1 + ci h, yn- 1 + y n+ 1 = ( 1 s b
∀ aK
j= 1
) y n + y n- 1 +
s
∀
s
biK n, i +
i= 1
∀
s
~
b iK n- 1, i,
i= 1
w here i∀ bi + i∀ b i = 1+ , ci = j∀ a ij, yn ~ y ( tn ), tn = nh, h > 0 is a step size and 0! = 1 = 1 = 1
~
It is the purpose of th is paper to investigate the stability properties of ( 2a) ~ ( 2c) w ith respect to th e fo l low ing delay differentia l equat ion: y ( t ) = Ly ( t) + M y ( t y ( t) = ! ( t)
0
t
0 , Re ( ∀) < 0 ;
t= 0 .
%
( 4)
That yie ld s a si m p le recurrence relation for the approx i m ations y n ( n = 1 , 2, # ) y n+ 1 = r( h ) y n + r ( h ) y n- 1 w here
[ 3]
2 The L Stability of the t w o step R unge K utta m ethods
Befo re actu ally dea ling w ith GPL stab ility o f tw o step Runge Kutta m ethods for de lay d ifferent ia l syste m. W e first study L stability o f t w o step Runge Ku tta m ethods for ODE. Apply ing ( 2a) ~ ( 2c) to the fo llow ing test equation : y ( t) = ∀ y ( t) y ( t) =
1
Introduction
W e consider the num erica l stab ility o f the in it ial va lu e problem fo r delay d ifferen tial equations : y ( t) = f ( t, y ( t), y ( t y ( t) = ( t) t! 0 , )) t 0 , ( 1)
% % %
T hus, if the polynom ial p ( z ) is a Schur po lynom ial for any h, Re( h ) < 0 , th e tw o step Runge Kutta m eth od ( 2a) ~ ( 2c) is A stab le for ODE. L et us study th e exponentia l so lu tion y ( t) = C e o f ( 4), where C is a constant to be deter m in ed and C & 0 , w e have # = ∀ . In particu lar , th e exponent ial so lution y ( t) = C e also satisfies y ( t+ h ) /y ( t ) = e . Consequently
蒋成香 : 两步 R unge K u tta方法求解延迟微分方程的 GPL 稳定性
345
In [ 1], Z. B artoszew sk i and Z. Jackiew icz stud ied the nu m er ic al stability o f the tw o step Runge Kutta m ethods for DDE s. Kuang and Cong et al . consid ered th e PL stab ility of the b lo ck - m ethods and the GPL stab ility of the Rosenbrock m ethods , respectiv e ly. In this paper , w e further stu dy th e asym ptot ic stab ility of tw o step Runge Kutta m ethods w ith respect to ( 3) . F irs, t the su ff ic ient condit io ns o f L stability for ord in ary differential equatio n( ODE ) are d iscussed . Then , th e m ain resu lt o f the paper is T heorem 3 . 1 in section 3 , wh ich states that the tw o step Runge Kutta m ethods ( 2a) ~ ( 2c) for DDE s is GPL stable if and only if the corresponding m ethods fo r ODE is L stable .
JIANG Cheng x iang
( Co llege of T ianH ua , Shangha iN o r m al U niversity , Shanghai 201815 paper deals w ith the num erical stab ility of the tw o step Runge Kutta m ethods for de lay d ifferent ia l equa tio ns( DDE s ) . T he conditio ns of L stability of tw o step Runge Kutta m ethods for ordin ary differentia l equation ( ODE ) are discussed , It is show n that tw o step Runge Ku tta m th ods is GPL stable for DDE s if and on ly if the correspondin g m ethods for ODE is L stable . K ey w ord s : de lay d ifferentia l equat io ns; t w o step Runge Kutta m ethod ; GPL stability CLC num ber : O 241 . 8 Docum en t code : A A rticle ID: 1000 5137( 2010) 04 0344 08
Received date : 2009 10 10 B iography: JI ANG Cheng x iang( 1979- ), m a le , lec turer , Co llege o f T ian H ua , Shanghai N or m a lU n iversity .
第 4期
T
)
t
0 ;
t! 0 ,
( 3)
w here y ( t) = (y 1 ( t), y 2 ( t ), #, yN ( t) ) , L and M are constant com plex N ∃N m atr ic es , > 0 , ! ( t ) denotes a g iv en vector valu e functio n and y ( t) is unknown for t> 0 . In [ 5], G. H u and T. M itisu i consid ered th e nu m er ic al stab ility of the Runge Kutta m ethods for DDE s.
~
% [ 13]
~
leads to that the po lynom ial p ( z) = z - r ( h ) z - r ( h ) is a Schur polynom ial if
2
%
~
r (h ) !
%
%
z - r( h )
fo r any z, z = 1 , fo r any z, z = 1 .
w here f, deno te g iv en functions , > 0 is a given constant and y ( t) is unknown for t > 0 . f satisfies conditions wh ich guarantee the ex istence of the unique so lu tion y to ( 1). For such condit ions , one can refer to [ 4] and [ 9]. M any m um erica l m ethods have been proposed for th e num erica l so lution o f the prob le m ( 1) . In th is pa per , w e are concerned w ith tw o step Runge Kutta m ethods( TSRK ) o f the for m K n, i = hf ( tn + ci h, yn +