D02 – Ordinary Differential Equations D02HBF – NAG Fortran Library Routine Document

合集下载

耦合常微分方程组

耦合常微分方程组耦合常微分方程组（coupled ordinary differential equations）是指由多个微分方程组成的系统，每个微分方程都与其他方程相关并相互影响。

这种方程组常出现在许多领域，如物理学、工程学和生物学等。

在本文中，我们将介绍耦合常微分方程组的基本概念、求解方法以及一些应用示例。

首先，让我们明确耦合常微分方程组的定义。

耦合常微分方程组由n个微分方程组成，每个方程都涉及n个未知函数和它们的导数。

这些方程之间存在相互关系，因此必须同时求解才能得到系统的完整解。

具体而言，耦合常微分方程组可以表示为:dx1/dt = f1(x1, x2, ..., xn)dx2/dt = f2(x1, x2, ..., xn)dxn/dt = fn(x1, x2, ..., xn)其中x1, x2, ..., xn是未知函数，f1, f2, ..., fn是给定的函数。

解决耦合常微分方程组的一种常见方法是数值求解。

数值求解可以将微分方程组转化为差分方程，然后使用数值方法逼近解。

其中一种常用的数值方法是欧拉法（Euler's method），它基于线性插值来逼近解。

欧拉法的基本思想是用导数值在某一点上的近似值来估计下一个点上的函数值。

通过迭代计算，可以逼近整个解曲线。

然而，欧拉法的精度有限，对于一些高阶的方程组可能不够准确。

因此，还有许多其他更高阶的数值方法可供选择，如龙格-库塔法（Runge-Kutta method）和Adams方法等。

另一种求解耦合常微分方程组的方法是使用解析解。

解析解是指用一种或多种已知的函数去表示解的形式。

通常，解析解只存在于特殊的情况下，对于大多数复杂的方程组而言，很难找到解析解。

然而，对于某些简单的方程组，可以通过分离变量、线性叠加或变量替换等方法，得到解析解。

解析解具有精确性和清晰的物理意义，因此在一些特定应用中很有用。

现在让我们来看一些实际的应用示例。

ode复习资料

ode复习资料ODE 复习资料ODE（Ordinary Differential Equations，常微分方程）是数学中的一个重要分支，研究的是描述自然界中变化现象的方程。

它在物理学、工程学、生物学等领域都有广泛的应用。

对于学习ODE的同学来说，复习资料是非常重要的辅助工具。

本文将介绍一些常见的ODE复习资料，并提供一些学习建议。

一、经典教材1. 《常微分方程教程》（作者：丁同仁）：这是一本经典的ODE教材，内容全面、系统，适合初学者。

书中深入浅出地介绍了ODE的基本概念、解法和应用，配有大量的例题和习题，有助于读者理解和巩固知识。

2. 《常微分方程教程》（作者：谢金星）：这本教材是国内外ODE教学的重要参考书之一。

作者结合自己多年的教学经验，将复杂的理论问题简化为易于理解的形式，让读者能够迅速掌握ODE的基本概念和解法。

3. 《Ordinary Differential Equations》（作者：Vladimir I. Arnold）：这是一本经典的英文教材，被誉为ODE领域的圣经。

书中内容丰富，涵盖了ODE的各个方面，包括基本理论、解的存在唯一性、稳定性等。

虽然是英文教材，但对于有一定数学基础的同学来说，是一本不可多得的学习资料。

二、在线课程1. Coursera：Coursera是一个知名的在线教育平台，提供了许多优质的ODE课程。

例如，由加州大学洛杉矶分校（UCLA）开设的《Ordinary Differential Equations for Engineers》课程，介绍了ODE的基本概念和解法，并通过实例讲解了ODE在工程领域的应用。

2. MIT OpenCourseWare：麻省理工学院的开放课程平台提供了许多免费的ODE课程资料。

例如，他们的《Differential Equations》课程包括了从ODE的基本概念到高级应用的内容，让学习者能够全面掌握ODE的知识。

三、学习建议1. 理论与实践结合：ODE是一门理论与实践相结合的学科，理论只有通过实际问题的解析才能更好地理解。

一阶常微分方程

Chapter 1 First-order ordinary differential equations (ODE)一階常微分方程1.1 基本概念()x f y =或()t f y =，y 是x 或t 的函數，y 是因變數(dependent variable )，x 或t 是自變數(independent variable )◎ 微分方程(differential equations)：一方程式包含有因變數y 關於自變數x 或t的導數(derivatives)y y ′ ,&或微分(differentials)dy 。

◎ 常微分方程(ordinary differential equations, ODE)：一微分方程包含有一個或數個因變數（通常為()x y ）關於僅有一個自變數x 的導數。

Ex. 222)2(2 ,09 ,cos y x y e y y x y y x y x +=′′+′′′′=+′′=′◎ 偏微分方程(partial differential equations, PDE)：一微分方程包含至少有一個因變數關於兩個以上自變數的部分導數。

Ex. 02222=∂∂+∂∂yux u◎ 微分方程的階數：在微分方程式中所出現最高階導數的階數。

◎ 線性微分方程：在微分方程式中所出現的因變數因變數因變數或其導數僅有一次式(first degree)而無二次以上的乘積（自變數可以有二次以上的乘積）。

Ex. x y y x y cos 24=+′+′′ 因變數：y ，自變數：x ，二階線性常微分方程 x y y y y cos 24=+′+′′ 因變數：y ，自變數：x ，二階非線性常微分方程 222)2(2 y x y e y y x x +=′′+′′′′ 因變數：y ，自變數：x ，三階非線性常微分方程□ 一階常微分方程(first-order ordinary differential equations)隱式形式(implicit form) 表示 0),,(=′y y x F (4)顯式形式(explicit form) 表示 ),(y x f y =′Ex. 隱式形式ODE 0423=−′−y y x ，當0≠x 時，可表示為顯式形式234y x y =′□ 解的概念(concept of solution)在某些開放間隔區間b x a <<，一函數)(x h y =是常微分方程常微分方程0),,(=′y y x F 的解，其函數)(x h 在此區間b x a <<是明確(defined)且可微分的(differentiable)，其)(x h 的曲線（或圖形）是被稱為解答曲線(solution curve)。

常微分方程的稳定解与不稳定解

常微分方程的稳定解与不稳定解常微分方程（Ordinary Differential Equations, ODE）是数学中重要的一门分支，研究函数的导数或微分在各种条件下的变化规律，广泛应用于物理、生物、工程等领域。

在解常微分方程的过程中，存在着两种重要的解：稳定解和不稳定解。

本文将对这两种解进行详细的介绍和分析。

1. 稳定解稳定解是指在一定条件下，系统的解向该解趋近，即当初始条件发生微小变化时，解会收敛到该解附近。

在常微分方程中，稳定解对应着系统的平衡点或稳定点，其解析形式通常为一组常数。

稳定解的性质可通过线性稳定性判据进行分析。

对于一阶常微分方程，即形如dy/dt = f(y)的方程，设y = c为方程的一个平衡解，则只需考虑f(c)的符号即可判断平衡解的稳定性：1.1 当f(c) < 0时，平衡解c是局部稳定解。

1.2 当f(c) > 0时，平衡解c是不稳定解。

例如，考虑一阶线性常微分方程dy/dt = -ky，其中k为正常数。

解析解为y = ce^(-kt)，其中c为常数。

当k > 0时，f(c) = -kc < 0，即平衡解y = 0是稳定解。

2. 不稳定解不稳定解指的是在一定条件下，系统的解远离该解，即当初始条件发生微小变化时，解会远离该解。

与稳定解相对应的，不稳定解对应着系统的不稳定点。

不稳定解的性质与稳定解相反，也可通过线性稳定性判据进行判断：2.1 当f(c) < 0时，平衡解c是不稳定解。

2.2 当f(c) > 0时，平衡解c是局部稳定解。

以二阶微分方程为例进行说明。

考虑二阶线性常微分方程d^2y/dt^2 + c1 * dy/dt + c2 * y = 0，其中c1和c2为常数。

该方程的解形式为y = Ae^(m1t) + Be^(m2t)，其中A和B为常数，m1和m2为方程的特征根。

根据特征根的性质，可判断解的稳定性：2.3 当特征根m1和m2的实部大于零时，平衡解是不稳定解。

常微分方程(王高雄)第三版 3.3

初值问题
dy f ( x, y ) , dx y ( x0 ) y0 (3.1) '
的解y ( x, x0 , y0 )都在区间 [a, b]上存在, 并且 ( x, x0 , y0 ) ( x, x0 , y0 ) , x [a, b] 则称初值问题(3.1) '的解y ( x, x0 , y0 )在点( x0 , y0 )
前提解存在唯一
y0 ( x0 , x, y )
证明在(3.1)满足y ( x0 ) y0的解存在区间内任取一值x1 ,
y1 ( x1 , x0 , y0 ), 则由解的唯一性知, (3.1)过点( x1 , y1 )与过点( x0 , y0 )的解是同一条积分曲线 , 即此解也可写成: y ( x, x1 , y1 ), 且显然有: y0 ( x0 , x1 , y1 ),
2 定理1 (解对初值的连续依赖性定理)
方程条件: I. f ( x , y ) 在G内连续且关于 y满足局部Lips.条件;
dy f ( x, y) , dx ( x, y) G R2 (1)
II. y ( x , x0 , y0 ) 是(1)满足( x0 , y0 ) G 的解,定义
C 时,有 S G G 覆盖定理,存在N,当G i i 1 对 0 ,记 y , S ), min , / 2 d (G
N
Ci
G
L max L1,, LN 则以为半径的圆,当其圆心从S的
G
左端点沿S 运动到右端点时,扫过的区域即为符合条件的要找区域D
0
义, 其中 a x0 b, 则对 0, ( , a, b) 0, 使当

常微分方程

若存在 ( x, c1 ,
, cn ) 的一个邻域，使得
, c1 ′ , c1 , c2 ′ , c2 , , cn ′ cn ( n 1) cn
≠0
( n 1) ( n 1) , , c1 c2
,
则称 y = ( x, c1 ,
, cn ) 含有n个相互独立的常数。
y 例： = c1 cos x + c2 sin x 是 y′′ + y = 0 的通解。因为 y′ = c1 sin x + c2 cos x 而
是
在 (∞, +∞)上的解。
2
y = tan(t )
例：xdx +
x = 1+ x
'
在 (
π π
, ) 上的解。 2 2
ydy = 0 有隐式解 x 2 + y 2 = C （ C 为任意常数）。
n 阶方程的通解：把含有 n 个相互独立的任意常数
c 称为 c1，c 2，， n 的解 y＝ x1，c1，，c n） n 阶方程的通解。（
耦合摆的动态演示
摆长减小的单摆
我们只研究这样一个方程:
θ( t ) 2 2 t 10 θ ( t ) + 2 θ( t ) =0 t 10 t 10 t
用Maple7编写的单摆模型的动态示意图
1.1.2 微分方程的基本概念
凡含有自变量、未知函数以及未知函数的导数（或微分）的方程称为微分方程。例如：
用maple 7解双摆的运动微分方程
2 2 θ1 ( t ) = 10 θ2 ( t ) 20 θ1 ( t ) t
2
2 2 θ2 ( t ) = 20 θ1 ( t ) 20 θ2 ( t ) t

neural ordinary differential equations生成式模型

neural ordinary differential equations生成式模型
神经常微分方程（Neural Ordinary Differential Equations，简称Neural ODE）是一种生成式模型，用于建模和生成时间序列
数据。

Neural ODE 的核心思想是使用常微分方程（Ordinary Differential Equations，简称 ODE）来建模数据的演化过程。

传统的生成模型通常基于离散时间步的迭代，而 Neural ODE
基于连续时间的积分形式，能够更好地捕捉数据的连续演化。

Neural ODE 的基本模型由两部分组成：ODE 函数和初始化条件。

ODE 函数定义了数据的演化规律，通常使用神经网络来
表示。

初始化条件为时间序列数据的初始状态。

生成时间序列数据的过程是通过解ODE函数得到的。

通过在
一定时间间隔内积分ODE函数，可以逐步推断出数据的演化
过程，并生成新的时间序列数据。

Neural ODE 在生成模型中的优势在于，它能够处理不规则、
不完整的时间序列数据，并且能够根据已有数据自适应地生成未来的数据。

此外，由于使用了ODE函数来建模数据的演化，生成的数据具有平滑性和连续性。

总之，Neural ODE 是一种生成式模型，通过使用常微分方程
来建模数据的演化过程，能够生成具有平滑性和连续性的时间序列数据，并适用于处理不规则和不完整的数据。

Differential Equations - Introduction

Part II Differential Equations - Introduction
Different first order differential equations
Standard form and differential form Standard form for a first-order differential equation in the unknown function y (t) is
Bernoulli’s equation
A Bernoulli differential equation is an equation of the form
Part II Differential Equations - Introduction
Different first order differential equations
PART II Differential Equations
Introduction
PART II Differential Equations - Introduction
Contents
Basic concepts Different first order differential equations
Separable equation Consider the differential form, if M(t, y) = A(t), and N(t, y) = B(y), the differential equation is separable.
Part II Differential Equations - Introduction
M(t, y)dt + N(t, y)dy = 0

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

D02–Ordinary Diﬀerential EquationsD02HBF–NAG Fortran Library Routine DocumentNote.Before using this routine,please read the Users’Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.1PurposeD02HBF solves the two-point boundary-value problem for a system of ordinary diﬀerential equations, using initial value techniques and Newton iteration;it generalises subroutine D02HAF to include the case where parameters other than boundary values are to be determined.2SpeciﬁcationSUBROUTINE D02HBF(P,N1,PE,E,N,SOLN,M1,FCN,BC,RANGE,W,1IW,IFAIL)INTEGER N1,N,M1,IW,IFAILreal P(N1),PE(N1),E(N),SOLN(N,M1),W(N,IW)EXTERNAL FCN,BC,RANGE3DescriptionD02HBF solves the two-point boundary-value problem by determining the unknown parametersp1,p2,...,p n1of the problem.These parameters may be,but need not be,boundary values;theymay include eigenvalue parameters in the coeﬃcients of the diﬀerential equations,length of the range of integration,etc.The notation and methods used are similar to those of D02HAF and the user is advised tostudy thisﬁrst.(The parameters p1,p2,...,p n1correspond precisely to the unknown boundary conditionsin D02HAF.)It is assumed that we have a system of nﬁrst-order ordinary diﬀerential equations of the form:dy idx=f i(x,y1,y2,...,y n),i=1,2,...,n,and that the derivatives fi are evaluated by a subroutine FCNsupplied by the user.The system,includingthe boundary conditions given by BC and the range of integration given by RANGE,involves the n1unknown parameters p1,p2,...,pn1which are to be determined,and for which initial estimates must besupplied.The number of unknown parameters n1must not exceed the number of equations n.If n1<n,we assume that(n−n1)equations of the system are not involved in the matching process.These are usually referred to as‘driving equations’;they are independent of the parameters and of the solutions of the other n1equations.In numbering the equations for the subroutine FCN,the driving equations must be putﬁrst.The estimated values of the parameters are corrected by a form of Newton iteration.The Newton correction on each iteration is calculated using a Jacobian matrix whose(i,j)th element depends on thederivative of the i th component of the solution,yi ,with respect to the j th parameter,pj.This matrix iscalculated by a simple numerical diﬀerentiation technique which requires n1evaluations of the diﬀerential system.If the parameter IFAIL is set appropriately,the routine automatically prints messages to inform the user of theﬂow of the calculation.These messages are discussed in detail in Section8.D02HBF is a simpliﬁed version of D02SAF which is described in detail in Gladwell[1].4References[1]Gladwell I(1979)The development of the boundary value codes in the ordinary diﬀerential equationschapter of the NAG Library Codes for Boundary Value Problems in Ordinary Diﬀerential Equations.Lecture Notes in Computer Science(ed B Childs,M Scott,J W Daniel,E Denman and P Nelson) 76Springer-Verlag5ParametersUsers are strongly recommended to read Section3and Section8in conjunction with this section.1:P(N1)—real array Input/OutputOn entry:an estimate for the i th parameter,pi ,for i=1,2,...,n1.On exit:the corrected value for the i th parameter,unless an error has occurred,when it contains the last calculated value of the parameter.2:N1—IN TEG ER Input On entry:the number of parameters,n1.Constraint:1≤N1≤N.3:PE(N1)—real array Input On entry:the elements of PE must be given small positive values.The element PE(i)is used(i)in the convergence test on the i th parameter in the Newton iteration,and(ii)in perturbing the i th parameter when approximating the derivatives of the components of the solution with respect to this parameter for use in the Newton iteration.The elements PE(i)should not be chosen too small.They should usually be several orders of magnitude larger than machine precision.Constraint:PE(i)>0.0,for i=1,2,...,N1.4:E(N)—real array Input On entry:the elements of E must be given positive values.The element E(i)is used in the bound on the local error in the i th component of the solution y i during integration.The elements E(i)should not be chosen too small.They should usually be several orders of magnitude larger than machine precision.Constraint:E(i)>0.0,for i=1,2,...,N.5:N—INTEGER Input On entry:the total number of diﬀerential equations,n.Constraint:N≥2.6:SOLN(N,M1)—real array Output On exit:the solution when M1>1(see below).7:M1—INTEGER Input On entry:a value which controls exit values as follows:M1=1Theﬁnal solution is not calculated;M1>1Theﬁnal values of the solution at interval(length of range)/(M1−1)are calculated and storedsequentially in the array SOLNstarting with the values of the solutions evaluated at theﬁrst end-point(see subroutine RANGE below)stored in theﬁrst column of SOLN.Constraint:M1≥1.FCNmust evaluate the function fi (i.e.,the derivative yi),for i=1,2,...,n.Its speciﬁcation is:SUBROUTINE FCN(X,Y,F,P)real X,Y(n),F(n),P(n1)where n and n1are the actual values of Nand N1in the call of D02HBF.1:X—real Input On entry:the value of the argument x.2:Y(n)—real array Input On entry:the value of the argument y i,for i=1,2,...,n.3:F(n)—real array Output On exit:the value of f i,for i=1,2,...,n.The f i may depend upon the parameters p j,forj=1,2,...,n1.If there are any driving equations(see Section3)then these must be numberedﬁrst in the ordering of the components of F in FCN.4:P(n1)—real array Input On entry:the current estimate of the parameter p i,for i=1,2,...,n1.FCNmust be declared as EXTERN AL in the(sub)program from which D02HBF is called.Parameters denoted as Input must not be changed by this procedure.9:BC—SUBROUTINE,supplied by the user.External Procedure BC must place in G1and G2the boundary conditions at a and b respectively(see RANGE below).Its speciﬁcation is:SUBROUTINE BC(G1,G2,P)real G1(n),G2(n),P(n1)where n and n1are the actual values of Nand N1in the call of D02HBF.1:G1(n)—real array Output On exit:the value of y i(a),(where this may be a known value or a function of the parametersp j ,for j=1,2,...,n1);i=1,2,...,n.2:G2(n)—real array Output On exit:the value of y i(b),for i=1,2,...,n,(where these may be known values or functionsof the parameters pj ,for j=1,2,...,n1).If n>n1,so that there are some driving equations,then theﬁrst n−n1values of G2need not be set since they are never used.3:P(n1)—real array Input On entry:an estimate of the parameter p i,for i=1,2,...,n1.BC must be declared as EXTERNAL in the(sub)program from which D02HBF is called.Parameters denoted as Input must not be changed by this procedure.RANGE must evaluate the boundary points a and b,each of which may depend on the parameters .The integrations in the shooting method are always from a to b.p1,p2,...,p n1Its speciﬁcation is:SUBROUTINE RANGE(A,B,P)real A,B,P(n1)where n1is the actual value of N1in the call of D02HBF.1:A—real Output On exit:one of the boundary points,a.2:B—real Output On exit:the second boundary point,b.N ote that B>A forces the direction of integration tobe that of increasing X.If A and B are interchanged the direction of integration is reversed.3:P(n1)—real array Input On entry:the current estimate of the i th parameter,p i,for i=1,2,...,n1.RANGE must be declared as EXTERNAL in the(sub)program from which D02HBF is called.Parameters denoted as Input must not be changed by this procedure.11:W(N,IW)—real array Output Used mainly for workspace.On exit:with IFAIL=2,3,4or5(see Section6),W(i,1),for i=1,2,...,n contains the solution at the point x when the error occurred.W(1,2)contains x.12:IW—INTEGER Input On entry:the second dimension of the array W as declared in the(sub)program from which D02HBF is called.Constraint:IW≥3N+14+max(11,N).13:IFAIL—INTEGER Input/Output For this routine,the normal use of IFAIL is extended to control the printing of error and warning messages as well as specifying hard or soft failure(see Chapter P01).Before entry,IFAIL must be set to a value with the decimal expansion cba,where each of the decimal digits c,b and a must have a value of0or1.a=0speciﬁes hard failure,otherwise soft failure;b=0suppresses error messages,otherwise error messages will be printed(see Section6);c=0suppresses warning messages,otherwise warning messages will be printed(see Section6).The recommended value for inexperienced users is110(i.e.,hard failure with all messages printed).Unless the routine detects an error(see Section6),IFAIL contains0on exit.6Error Indicators and WarningsFor each error,an explanatory error message is output on the current error message unit(as deﬁned by X04AAF),unless suppressed by the value of IFAIL on entry.Errors detected by the routine:IFAIL=1One or more of the parameters N,N1,M1,IW,E or PE is incorrectly set.IFAIL=2The step length for the integration became too short whilst calculating the residual(see Section8).IFAIL=3No initial step length could be chosen for the integration whilst calculating the residual.Note:IFAIL=2or3can occur due to choosing too small a value for E or due to choosing the wrong direction of integration.Try varying E and interchanging a and b.These error exits can also occur for very poor initial choices of the parameters in the array P and,in extreme cases, because this routine cannot be used to solve the problem posed.IFAIL=4As for IFAIL=2but the error occurred when calculating the Jacobian.IFAIL=5As for IFAIL=3but the error occurred when calculating the Jacobian.IFAIL=6The calculated Jacobian has an insigniﬁcant column.This can occur because a parameter p i is incorrectly entered when posing the problem.Note:IFAIL=4,5or6usually indicate a badly scaled problem.The user may vary the size of PE.Otherwise the use of the more general D02SAF which aﬀords more control over the calculations is advised.IFAIL=7The linear algebra routine used(F02WEF)has failed.This error exit should not occur and can be avoided by changing the initial estimates p.iIFAIL=8The Newton iteration has failed to converge.This can indicate a poor initial choice of parameters p i or a very diﬃcult problem.Consider varying the elements PE(i)if the residuals are small in the monitoring output.If the residuals are large,try varying the initial parameters p.iIFAIL=9,10,11,12or13Indicate that a serious error has occurred in D02SAZ,D02SAW,D02SAX,D02SAU or D02SAV respectively.Check all array subscripts and subroutine parameter lists in the call to D02HBF.Seek expert help.7AccuracyIf the process converges,the accuracy to which the unknown parameters are determined is usually close to that speciﬁed by the user;and the solution,if requested,may be determined to a required accuracy by varying the parameter E.8Further CommentsThe time taken by the routine depends on the complexity of the system,and on the number of iterations required.In practice,integration of the diﬀerential equations is by far the most costly process involved. Wherever they occur in the routine,the error parameters contained in the arrays E and PE are used in ‘mixed’form;that is E(i)always occurs in expressions of the formE(i)×(1+|y i|)and PE(i)always occurs in expressions of the formPE(i)×(1+|p i|)Though not ideal for every application,it is expected that this mixture of absolute and relative error testing will be adequate for most purposes.The user may determine a suitable direction of integration a to b and suitable values for E(i)by integrations with D02PCF.The best direction of integration is usually the direction of decreasing solutions.The user is strongly recommended to set IFAIL to obtain self-explanatory error messages, and also monitoring information about the course of the computation.The user may select the channel numbers on which this output is to appear by calls of X04AAF(for error messages)or X04ABF(for monitoring information)–see Section9for an example.Otherwise the default channel numbers will be used,as speciﬁed in the implementation document.The monitoring information produced at each iteration includes the current parameter values,the residuals and two norms:a basic norm and a current norm.At each iteration the aim is toﬁnd parameter values which make the current norm less than the basic norm.Both these norms should tend to zero as should the residuals.(They would all be zero if the exact parameters were used as input.)For more details,in particular about the other monitoring information printed,the user is advised to consult the speciﬁcation of D02SAF and,especially,the description of the parameter MONIT there.The computing time for integrating the diﬀerential equations can sometimes depend critically on thequality of the initial estimates for the parameters pi .If it seems that too much computing time isrequired and,in particular,if the values of the residuals printed by the monitoring routine are much larger than the expected values of the solution at b then the coding of the subroutines FCN,BC and RANGE should be checked for errors.If no errors can be found,an independent attempt should be made to improve the initial estimates for p i.The subroutine can be used to solve a very wide range of problems,for example:(a)eigenvalue problems,including problems where the eigenvalue occurs in the boundary conditions;(b)problems where the diﬀerential equations depend on some parameters which are to be determinedso as to satisfy certain boundary conditions(see example(ii)in Section9);(c)problems where one of the end-points of the range of integration is to be determined as the pointwhere a variable y i takes a particular value(see example(ii)in Section9);(d)singular problems and problems on inﬁnite ranges of integration where the values of the solution ata orb or both are determined by a power series or an asymptotic expansion(or a more complicatedexpression)and where some of the coeﬃcients in the expression are to be determined(see example(i)in Section9);and(e)diﬀerential equations with certain terms deﬁned by other independent(driving)diﬀerentialequations.9ExampleFor this routine two examples are presented,in Section9.1and Section9.2.In the example programs distributed to sites,there is a single example program for D02HBF,with a main program: *D02HBF Example Program Text*Mark14Revised.NAG Copyright1989.*..Parameters..INTEGER NOUTPARAMETER(NOUT=6)*..External Subroutines..EXTERNAL EX1,EX2*..Executable Statements..WRITE(NOUT,*)’D02HBF Example Program Results’CALL EX1CALL EX2STOPENDThe code to solve the two example problems is given in the subroutines EX1and EX2,in Section9.1.1 and Section9.2.1respectively.9.1Example1Toﬁnd the solution of the diﬀerential equationy =(y3−y )/2xon the range0≤x≤16,with boundary conditions y(0)=0.1and y(16)=1/6.We cannot use the diﬀerential equation at x=0because it is singular,so we take a truncated power series expansiony(x)=1/10+p1×√x/10+x/100near the origin where p1is one of the parameters to be determined.We choose the interval as[0.1,16] and setting p2=y (16),we can determine all the boundary conditions.We take X1=16.We write y =Y(1),y =Y(2),and estimate PARAM(1)=0.2,PARAM(2)=0.0.Note the call to X04ABF before the call to D02HBF.9.1.1Program TextNote.The listing of the example program presented below uses bold italicised terms to denote precision-dependent details. Please read the Users’Note for your implementation to check the interpretation of these terms.As explained in the Essential Introduction to this manual,the results produced may not be identical for all implementations.*SUBROUTINE EX1*..Parameters..INTEGER NOUTPARAMETER(NOUT=6)INTEGER N,N1,IW,M1PARAMETER(N=2,N1=2,IW=3*N+14+11,M1=6)*..Local Scalars..real H,X,X1INTEGER I,IFAIL,J*..Local Arrays..real C(N,M1),ERROR(N),PARAM(N1),PARERR(N1),W(N,IW) *..External Subroutines..EXTERNAL AUX1,BCAUX1,D02HBF,RNAUX1,X04ABF*..Intrinsic Functions..INTRINSIC real*..Executable Statements..WRITE(NOUT,*)WRITE(NOUT,*)WRITE(NOUT,*)’Case1’CALL X04ABF(1,NOUT)PARAM(1)=0.2e0PARAM(2)=0.0e0PARERR(1)=1.0e-5PARERR(2)=1.0e-3ERROR(1)=1.0e-4ERROR(2)=1.0e-4**Set IFAIL to111to obtain monitoring information*IFAIL=11*CALL D02HBF(PARAM,N1,PARERR,ERROR,N,C,M1,AUX1,BCAUX1,RNAUX1,W,IW,+IFAIL)*WRITE(NOUT,*)IF(IFAIL.NE.0)THENWRITE(NOUT,99999)’IFAIL=’,IFAILIF(IFAIL.LE.5.AND.IFAIL.NE.1)THENWRITE(NOUT,*)WRITE(NOUT,99996)’W(1,2)=’,W(1,2),’W(.,1)=’, +(W(I,1),I=1,N)END IFELSEWRITE(NOUT,*)’Final parameters’WRITE(NOUT,99998)(PARAM(I),I=1,N1)WRITE(NOUT,*)WRITE(NOUT,*)’Final solution’WRITE(NOUT,*)’X-value Components of solution’CALL RNAUX1(X,X1,PARAM)H=(X1-X)/real(M1-1)DO20I=1,M1WRITE(NOUT,99997)X+(I-1)*H,(C(J,I),J=1,N) 20CONTINUEEND IFRETURN*99999FORMAT(1X,A,I3)99998FORMAT(1X,1P,3e15.3)99997FORMAT(1X,F7.2,2F13.4)99996FORMAT(1X,A,F9.4,A,10e10.3)END*SUBROUTINE AUX1(X,Y,F,PARAM)*..Parameters..INTEGER NPARAMETER(N=2)*..Scalar Arguments..real X*..Array Arguments..real F(N),PARAM(N),Y(N)*..Executable Statements..F(1)=Y(2)F(2)=(Y(1)**3-Y(2))/(2.0e0*X)RETURNEND*SUBROUTINE RNAUX1(X,X1,PARAM)*..Scalar Arguments..real X,X1*..Array Arguments..real PARAM(2)*..Executable Statements..X=0.1e0X1=16.0e0RETURNEND*SUBROUTINE BCAUX1(G,G1,PARAM)*..Parameters..INTEGER NPARAMETER(N=2)*..Array Arguments..real G(N),G1(N),PARAM(N)*..Local Scalars..real Z*..Intrinsic Functions..INTRINSIC SQRT*..Executable Statements..Z=0.1e0G(1)=0.1e0+PARAM(1)*SQRT(Z)*0.1e0+0.01e0*ZG(2)=PARAM(1)*0.05e0/SQRT(Z)+0.01e0G1(1)=1.0e0/6.0e0G1(2)=PARAM(2)RETURNEND9.1.2Program DataNone.9.1.3Program ResultsD02HBF Example Program ResultsCase1Final parameters4.629E-02 3.494E-03Final solutionX-value Components of solution0.100.10250.01733.280.12170.00426.460.13380.00369.640.14490.003412.820.15570.003416.000.16670.00359.2Example2Toﬁnd the gravitational constant p1and the range p2over which a projectile must beﬁred to hit the target with a given velocity.The diﬀerential equations arey =tanφv =−(p1sinφ+0.00002v2)v cosφφ =−p1 v2on the range0<x<p2,with boundary conditionsy=0,v=500,φ=0.5at x=0,y=0,v=450,φ=p3at x=p2.We write y=Y(1),v=Y(2),φ=Y(3).We estimate p1=PARAM(1)=32,p2=PARAM(2)=6000and p3=PARAM(3)=0.54(though this last estimate is not important).9.2.1Program TextNote.The listing of the example program presented below uses bold italicised terms to denote precision-dependent details. Please read the Users’Note for your implementation to check the interpretation of these terms.As explained in the Essential Introduction to this manual,the results produced may not be identical for all implementations.*SUBROUTINE EX2*..Parameters..INTEGER NOUTPARAMETER(NOUT=6)INTEGER N,N1,IW,M1PARAMETER(N=3,N1=3,IW=3*N+14+11,M1=6)*..Local Scalars..real H,X,X1INTEGER I,IFAIL,J*..Local Arrays..real C(N,M1),ERROR(N),PARAM(N1),PARERR(N1),W(N,IW) *..External Subroutines..EXTERNAL AUX2,BCAUX2,D02HBF,RNAUX2,X04ABF*..Intrinsic Functions..INTRINSIC real*..Executable Statements..WRITE(NOUT,*)WRITE(NOUT,*)WRITE(NOUT,*)’Case2’CALL X04ABF(1,NOUT)PARAM(1)=32.0e0PARAM(2)=6000.0e0PARAM(3)=0.54e0PARERR(1)=1.0e-5PARERR(2)=1.0e-4PARERR(3)=1.0e-4ERROR(1)=1.0e-2ERROR(2)=1.0e-2ERROR(3)=1.0e-2**Set IFAIL to111to obtain monitoring information*IFAIL=11*CALL D02HBF(PARAM,N1,PARERR,ERROR,N,C,M1,AUX2,BCAUX2,RNAUX2,W,IW,+IFAIL)*WRITE(NOUT,*)IF(IFAIL.NE.0)THENWRITE(NOUT,99999)’IFAIL=’,IFAILIF(IFAIL.LE.5.AND.IFAIL.NE.1)THENWRITE(NOUT,*)WRITE(NOUT,99996)’W(1,2)=’,W(1,2),’W(.,1)=’,+(W(I,1),I=1,N)END IFELSEWRITE(NOUT,*)’Final parameters’WRITE(NOUT,99998)(PARAM(I),I=1,N1)WRITE(NOUT,*)WRITE(NOUT,*)’Final solution’WRITE(NOUT,*)’X-value Components of solution’CALL RNAUX2(X,X1,PARAM)H=(X1-X)/real(M1-1)DO20I=1,M1WRITE(NOUT,99997)X+(I-1)*H,(C(J,I),J=1,N)20CONTINUEEND IFRETURN*99999FORMAT(1X,A,I3)99998FORMAT(1X,1P,3e15.3)99997FORMAT(1X,F7.0,2F13.1,F13.3)99996FORMAT(1X,A,F9.4,A,10e10.3)END*SUBROUTINE AUX2(X,Y,F,PARAM)*..Parameters..INTEGER NPARAMETER(N=3)*..Scalar Arguments..real X*..Array Arguments..real F(N),PARAM(N),Y(N)*..Intrinsic Functions..INTRINSIC COS,TAN*..Executable Statements..F(1)=TAN(Y(3))F(2)=-PARAM(1)*TAN(Y(3))/Y(2)-0.00002e0*Y(2)/COS(Y(3))F(3)=-PARAM(1)/Y(2)**2RETURNEND*SUBROUTINE RNAUX2(X,X1,PARAM)*..Parameters..INTEGER NPARAMETER(N=3)*..Scalar Arguments..real X,X1*..Array Arguments..real PARAM(N)*..Executable Statements..X=0.0e0X1=PARAM(2)RETURNEND*SUBROUTINE BCAUX2(G,G1,PARAM)*..Parameters..INTEGER NPARAMETER(N=3)*..Array Arguments..real G(N),G1(N),PARAM(N)*..Executable Statements..G(1)=0.0e0G(2)=500.0e0G(3)=0.5e0G1(1)=0.0e0G1(2)=450.0e0G1(3)=PARAM(3)RETURNEND[NP3390/19/pdf]D02HBF.119.2.2Program DataNone.9.2.3Program ResultsCase2Final parameters3.239E+01 5.962E+03-5.353E-01Final solutionX-value Components of solution0.0.0500.00.5001192.529.6451.60.3282385.807.2420.30.1233577.820.4409.4-0.1034769.556.1420.0-0.3305962.0.0450.0-0.535D02HBF.12(last)[NP3390/19/pdf]。