Hybrid Symbolic-Numeric Methods Applied to Definite Integrals and ODEsK.O. Geddes
Symbolic Computation Group
Department of Computer Science
University of Waterloo
Waterloo ON N2L 3G1
CANADA
http://daisy.uwaterloo.ca/~kogeddes<Text-field layout="Heading 1" style="Heading 1">Abstract</Text-field>A scientific computation environment such as Maple supports both symbolic and numeric mathematical computation. By exploiting the paradigm of hybrid symbolic-numeric computation, an enhanced level of computational power can be achieved in various problem areas. The strategy is to apply an appropriate combination of symbolic mathematical analysis and numerical computation.Hybrid symbolic-numeric methods are applied in Maple for the evaluation of definite integrals in the presence of singularities. For the case of solving ordinary differential equations, we have been investigating techniques to generate numerical update formulas for initial-value problems by applying a symbolic calculus of vector divided differences suggested by W. Kahan.<Text-field layout="Heading 1" style="Heading 1">Hybrid Symbolic-Numeric Integration</Text-field><Text-field layout="Heading 2" style="Heading 2">Example 1</Text-field>intf := Int(exp(-x^3)/(x^2+1), x=0..1);The integral does not evaluate in symbolic mode.value(intf);Evaluate to a numerical value via numerical integration.evalf(intf);Or to higher precisionevalf(intf, 25);<Text-field layout="Heading 2" style="Heading 2">Example 2</Text-field>g := sin(x)*ln(x)*exp(-x^3):intg := int(g, x=0..infinity);Maple is not able to express the definite integral in closed form.Evaluating numerically (to 25 digits) yields the following result.evalf(intg, 25);In addition to dealing with the infinite interval of integration, note that the integrand has a logarithmic singularity at x=0 as the following series expansion shows.series(g, x=0);The hybrid symbolic-numeric techniques being used by Maple allow reasonably efficient calculation even when very high precision is requested. The computation of the above numerical results employs several of the techniques discussed in the following sections. For more detailed presentations of this approach see [Geddes86], [Geddes92].<Text-field layout="Heading 2" style="Heading 2">The Case of an Analytic Integrand on a Finite Interval</Text-field>easy case: finite interval of integration, and the integrand is an analytic function with no singularities on the closed intervalmost traditional numerical integration methods could serve for this purposehowever we want methods which are very efficient, including computations withvery high precision<Text-field layout="Heading 3" style="Heading 3">Clenshaw-Curtis quadrature</Text-field>preferred default method in Maple, named NiMlJ2NjcXVhZEc=based on computing Chebyshev series expansion of integrand on specified interval and then integrating term-by-termNiMlJ2NjcXVhZEc= at its best when the integrand is analytic in a sufficiently large region of the complex plane surrounding the interval of integration, in which case the Chebyshev series expansion converges quickly<Text-field layout="Heading 3" style="Heading 3">Current Implementation Technique in Maple </Text-field>NiMlJ2NjcXVhZEc= is programmed to detect when convergence is too slowslow convergence implies that there is a nearby singularitywhen this occurs, return to the control routine which will invoke the singularity-handling routinesif attempts to remove the singularity are unsuccessful, invoke NiMlKHF1YWRleHBH, an adaptive double-exponential quadrature methodNiMlKHF1YWRleHBH is less sensitive to the presence of nearby singularities, or even singularities at the endpoints of the interval<Text-field layout="Heading 2" style="Heading 2">Generalized Series Expansions</Text-field>main tool for analyzing integrand singularities: generalized series expansionsTaylor series, or Laurent series with finite principal part, if it existsmore generally, not only powers of NiMlInhH but also terms involving singular functions such asNiMtJSNsbkc2IyUieEc=, NiMtJSRleHBHNiMsJComIiIiRiglInhHISIiRio=, NiMtJSRleHBHNiMsJComIiIiRigqJCUieEciIiMhIiJGLA==, etc.mathematical validity of such series relies on the fact that the non-polynomial functions appearing in the expansion can be bounded by powers of NiMlInhHfor example, for all NiMlInhH > 0NiMsJComIiIiRiUlInhHISIiRic=
< NiMtJSNsbkc2IyUieEc= < NiMlInhHNiMiIiE=
< NiMtJSRleHBHNiMsJComIiIiRiglInhHISIiRio= < NiMiIiI=NiMiIiE=
< NiMtJSRleHBHNiMsJComIiIiRigqJCUieEciIiMhIiJGLA== < 1References: Very powerful algorithms for computing generalized series expansions have been developed. Some early work was presented in [GeddesGonnet89]. For more recent algorithmic developments see [Salvy91], [Salvy92], [Richardson96].<Text-field layout="Heading 2" style="Heading 2">Example 3: Generalized series expansions</Text-field>series(sqrt(sin(x)), x=0, 5);series( ln(1-cos(2*x)), x=0, 8 );`evalf/int/genseries`( 1/(1-x*exp(-cos(x)/x)), x, 6):convert(%,polynom);Note: The latter example is a case which is not currently handled by Maple's series command but rather it is generated by a specialized routine.<Text-field layout="Heading 2" style="Heading 2">Transforming to an Analytic Integrand</Text-field>if integrand has a singularity at a point of the interval of integration, look for a transformation which will remove the singularitytwo types of transformations:
- subtracting off the singularity
- algebraic transformation of variables<Text-field layout="Heading 2" style="Heading 2">Example 4: Subtracting off the singularity</Text-field>h := ln(1 - cos(2*x));Int(h, x=0..1);Note that the graph of this integrand goes to NiMsJCUpaW5maW5pdHlHISIi at NiMvJSJ4RyIiIQ==. However, this is an integrable singularity - it behaves like NiMtJSNsbkc2IyUieEc= near NiMvJSJ4RyIiIQ==.plot(h, x=0..1);The generalized series expansion of h at x = 0 takes the following form.series(h, x=0, 8);The non-regular part isq := 2*ln(x);The new expressionnewh := h - q;is analytic on the interval [0, 1].Thus it can be integrated easily by ccquad.r1 := evalf(Int(newh, x=0..1));(computing to 10 digits of accuracy).Integrating q is easy because it has the indefinite integralint(q, x);and its definite integral can therefore be computed symbolically.r2 := int(q, x=0..1);Finally, summing the two values, we obtain the value for the original definite integration problem.Int(h, x=0..1) = r1 + r2;<Text-field layout="Heading 2" style="Heading 2">Example 5: Algebraic transformation of variables</Text-field>F := sqrt(sin(x));Int(F, x=0..2);The generalized series expansion of F at x = 0 is of the formseries(F, x=0, 5);Applying the change of variables t = x^(1/2) yieldsr3 := Int(eval(F,x=t^2) * diff(t^2,t), t = 0..sqrt(2));The new integrand is analytic and therefore it can be integrated easily by ccquad.evalf(r3);<Text-field layout="Heading 2" style="Heading 2">Direct Integration of a Generalized Series</Text-field>integrand has a singularity, suppose it is at the left end-point NiMvJSJ4RyIiIQ==no transformation found to remove the singularitytherefore, treat the generalized series expansion as an approximation of the integrand near the singularity and directly integrate this series term-by-term over its interval of accuracyremainder of the interval can then be handled by the numerical integration methodit is important to move as far as possible away from the singularity before resuming with the numerical method, because the efficiency of numerical methods is affected by the nearness of a singularitysuccess of this technique relies on the power of the symbolic integrator to handle many of the singular functions which may arisefor example, term-by-term integration of the generalized series might involve integrals of the form:int( ln(x)*x^2, x );int( exp(-1/x)*x^2, x );int( exp(-1/x^2)*x^2, x );<Text-field layout="Heading 2" style="Heading 2">Example 6</Text-field>G := exp(v - v^2/2) / (1 + 1/2*exp(v));Int(G, v=0..infinity);First, the interval is split into NiM7IiIhIiIi and NiM7IiIiJSlpbmZpbml0eUc= .For the finite interval, ccquad is applied directly.r01 := evalf(Int(G, v=0..1));For the infinite interval, the change of variables NiMvJSJ2RyomIiIiRiYlInhHISIi transforms the problem into the new integration problemr1inf := Int(-eval(G,v=1/x) * diff(1/x,x), x = 0..1);Let the integrand appearing here be named g. The first few terms of the generalized series expansion of g are as follows.g := op(1,r1inf);s := `evalf/int/genseries`(g, x, 6);Maple's symbolic integrator determines that the first two terms of s have the following indefinite integrals.terms := [seq(simplify(op(i,s),symbolic), i=1..6)]:int(terms[1],x);int(terms[2],x);Similarly, the integral of each term is computed.If we compute the definite integral over [0, 0.25], say, of the successive terms of the generalized series expansion of g, we find that the successive values are as follows.Digits := 20:seq( int(terms[i], x=0..0.25), i = 1..6 ):evalf( [ % ] );Clearly, the series representation is rapidly converging on this interval.Summing up these values yields the following result for the integralof g over [0, 0.25].r1 := convert( %, `+` );Digits := 10:Look at a plot of the integrand NiNJImdHNiI= and its series approximation NiNJInNHNiI= over the interval NiM3JCIiISIiIg== .Agreement between the two can be seen on NiM3JCIiIS1JJkZsb2F0R0kqcHJvdGVjdGVkR0YnNiQiI0QhIiM= and perhaps even further to about NiMtSSZGbG9hdEdJKnByb3RlY3RlZEdGJTYkIiImISIi (to a lower degree of accuracy).s := convert(s, 'polynom'):plot({g,s}, x=0..1);For the remaining interval [0.25, 1], ordinary numerical integration methods encounter no difficulties because there are no nearby singularities.r2 := evalf(Int(g, x=0.25..1));Sum these two values to get the integral r1inf.r1inf;r1inf := r1 + r2;Finally, we have obtained the answer for the original integration problem.Int(G, v=0..infinity) = r01 + r1inf;<Text-field layout="Heading 1" style="Heading 1">Exploiting the Hybrid Symbolic-Numeric Paradigm for ODEs</Text-field>A natural question to ask is: In what manner could the hybrid symbolic-numeric paradigm, which was successfully applied to the case of numerical definite integration, be applied to the numerical solution of ODEs? Investigations on this question are in progress.The paper [Geddes99] presents an investigation of a calculus of divided differences of vector functions, proposed by W. Kahan [Kahan93]. For a given ODE initial-value problem, it is possible to generate numerical update formulas for the system of ODEs by applying the divided difference calculus. The general concept is to exploit the power of symbolic computation to generate a numerical update formula which is particularly tailored to the given problem. In some cases, the resulting numerical computation phase can achieve enhanced accuracy and efficiency compared with traditional numerical methods.<Text-field layout="Heading 1" style="Heading 1">A Class of Numerical ODE Solvers</Text-field><Text-field layout="Heading 2" style="Heading 2">The Initial Value Problem</Text-field>Consider an initial-value problem (IVP)NiNDJC8tJSVkaWZmRzYkLSUieUc2IyUidEdGKy0lImZHNiNGKS8tJSh+fn5+fn55RzYjJkYrNiMiIiEmRilGNA==where NiMtJSJ5RzYjJSJ0Rw== is a vector and NiMtJSJmRzYjJSJ5Rw== is a vector function.Typically the right hand side might be written as NiMtJSJmRzYkJSJ0RyUieUc= indicating explicit dependence on both NiMlInRH and NiMlInlH. However, it is always possible to convert such a problem into the form of an autonomous initial value problem (AIVP), and we assume this form here.<Text-field layout="Heading 2" style="Heading 2">Numerical Update Formulas</Text-field>To obtain a numerical solution, one determines an update formulaNiMvJiUieUc2IywmJSJuRyIiIkYpRiktJSJRRzYkJSJoRyZGJTYjRig=to advance the (approximate) solution from NiMtJSJ5RzYjJiUidEc2IyUibkc= to NiMvLSUieUc2IyYlInRHNiMsJiUibkciIiJGLEYsLUYlNiMsJiZGKDYjRitGLCUiaEdGLA== for some chosen step size NiMlImhH .<Text-field layout="Heading 3" style="Heading 3">Euler's Method</Text-field>NiMvLSUiUUc2JCUiaEclImdHLCZGKCIiIiomRidGKi0lImZHNiNGKEYqRio=is a method with first-order convergence. It is derived from a Taylor series expansion through first-order terms. This formula suffers from two disadvantages as a numerical solver: low-order convergence; and it is an explicit formula (which has implications for numerical stability).<Text-field layout="Heading 3" style="Heading 3">Euler's Backward Formula</Text-field>Solve NiMvJSJRRywmJSJnRyIiIiomJSJoR0YnLSUiZkc2I0YkRidGJw== for NiMvJSJRRy1GJDYkJSJoRyUiZ0c= .This formula also has first-order convergence but it is an implicit formula. Implicit formulas have important advantages for numerical stability, but they imply additional computational steps each iteration.<Text-field layout="Heading 3" style="Heading 3">Second Order Formula based on the Jacobian</Text-field>Corresponding to Euler's Method which was based on a Taylor series expansion through first-order terms, we get the following update formula from a Taylor series expansion through second-order terms.NiMvLSUiUUc2JCUiaEclImdHLChGKCIiIiomRidGKi0lImZHNiNGKEYqRioqKC0lIUc2IyomRiciIiNGNCEiIkYqLSYlIkpHNiNGLUYuRipGLEYqRio=where NiMmJSJKRzYjJSJmRw== denotes the Jacobian matrix of NiMlImZH . This formula has second-order convergence but it is an explicit formula.Not discussed here is the question of obtaining higher-order formulas (promotion of second-order formulas to fourth-order, or higher order). See [Kahan97b] for an interesting approach.The approach we take here to turn the above second-order explicit formula into an implicit formula is to replace the Jacobian matrix NiMtJiUiSkc2IyUiZkc2IyUiZ0c= by a divided difference matrix NiMtJSREREZHNiQlImdHJSJRRw== which approximates NiMtJiUiSkc2IyUiZkc2IyUiZ0c= for nearby vectors NiMlImdH and NiMlIlFH . This approach has some advantages in its ability to generate a numerical update formula tailored to the given initial-value problem. It will be seen that a computer algebra system is indisensible for performing the algebra required by this technique.<Text-field layout="Heading 1" style="Heading 1">Divided Differences of Vector Functions</Text-field>We wish to define a matrix which will approximate the Jacobian matrix NiMtJiUiSkc2IyUiZkc2IyUiZ0c= . Just as a scalar derivative at a specified point can be approximated by a scalar divided difference based on evaluations at two "nearby" points, we define a matrixNiMtJSREREZHNiQlIlhHJSJaRw==(a "divided difference function") which satisfies the propertyNiMvLCYtJSh+fn5+fn5mRzYjJSJYRyIiIi0lImZHNiMlIlpHISIiKiYtJSREREZHNiRGKEYtRiksJkYoRilGLUYuRik= .Note that NiMlIlhH and NiMlIlpH are vectors.The NiMlJERERkc= matrix is not uniquely determined by the above property unless the dimension of the vector space is 1. In the 1-dimensional case,NiMvLSUkRERGRzYkJSJYRyUiWkcqJiwmLSUiZkc2I0YnIiIiLUYsNiNGKCEiIkYuLCZGJ0YuRihGMUYx and it equals NiMtJSVkaWZmRzYkLSUiZkc2IyUiWkdGKQ== when NiMvJSJYRyUiWkc= .In the n-dimensional case, one possible definition of the matrix NiMtJSREREZHNiQlIlhHJSJaRw== isThis is not the preferred definition because the integration can introduce new transcendental functions (e.g. logarithms) and, in general, it yields formulas which are more complicated than necessary. As pointed out by Kahan, it is preferable to consider a dog-leg path of integration consisting of segments parallel to coordinate axes. The resulting matrix NiMtJSREREZHNiQlIlhHJSJaRw== is composed of entries which are not much more complicated than the original function NiMlImZH ; the entries are simply scalar divided differences involving components of the vectors NiMlIlhH and NiMlIlpH .<Text-field layout="Heading 2" style="Heading 2">Example 1a: Straight-line Integration</Text-field>Consider the system of three equations NiMvLSUjeSdHNiMlInRHLSUiZkc2IyUieUc= where restart;with(LinearAlgebra):y := <y1, y2, y3>:f := v -> <1, v[3], cos(v[1])*v[3]+(1-cos(v[1]))*v[2]>:'f(y)' = f(y);The Jacobian matrix isJ := linalg[jacobian](f(y), y);The straight-line integral definition for NiMtJSREREZHNiQlIlhHJSJaRw== would yield a matrix whose first two rows are trivial, but the third row of the matrix consists of nontrivial expressions. The most complicated expression is the (3,1) entry and the result of straight-line integration of this entry is as follows.X := <X1, X2, X3>: Z := <Z1, Z2, Z3>:J31 := unapply(J[3,1], y1,y2,y3);int(J31(X1+theta*(Z1-X1), X2+theta*(Z2-X2), X3+theta*(Z3-X3)), theta=0..1);In contrast, the corresponding expression computed in Example 1b below is much simpler.<Text-field layout="Heading 2" style="Heading 2">Example 1b: Integration via a Dog-leg Path</Text-field>In the three-dimensional case, we wish to integrate from a point NiMvJSJYRy0lJ3ZlY3Rvckc2IzclJSNYMUclI1gyRyUjWDNH to a point NiMvJSJaRy0lJ3ZlY3Rvckc2IzclJSNaMUclI1oyRyUjWjNH. We choose to follow a dog-leg path consisting of segments parallel to the coordinate axes so that integration along any particular segment will be single-variable integration. As we prove in the following section, the NiMlImpH-NiMlI3RoRw== column of the matrix NiMtJSREREZHNiQlIlhHJSJaRw== can be defined by scalar divided differences in the NiMlImpH-NiMlI3RoRw== variable.g1 := v -> f(v)[1]:
g2 := v -> f(v)[2]:
g3 := v -> f(v)[3]:DDF := Matrix(3,3):DDF[1,1] := ( g1(<Z[1],X[2],X[3]>) - g1(X) ) / (Z[1] - X[1]):
DDF[2,1] := ( g2(<Z[1],X[2],X[3]>) - g2(X) ) / (Z[1] - X[1]):
DDF[3,1] := ( g3(<Z[1],X[2],X[3]>) - g3(X) ) / (Z[1] - X[1]):
ZXX := <Z[1],X[2],X[3]>:
DDF[1,2] := ( g1(<Z[1],Z[2],X[3]>) - g1(ZXX) ) / (Z[2] - X[2]):
DDF[2,2] := ( g2(<Z[1],Z[2],X[3]>) - g2(ZXX) ) / (Z[2] - X[2]):
DDF[3,2] := ( g3(<Z[1],Z[2],X[3]>) - g3(ZXX) ) / (Z[2] - X[2]):
DDF[1,3] := ( g1(Z) - g1(<Z[1],Z[2],X[3]>) ) / (Z[3] - X[3]):
DDF[2,3] := ( g2(Z) - g2(<Z[1],Z[2],X[3]>) ) / (Z[3] - X[3]):
DDF[3,3] := ( g3(Z) - g3(<Z[1],Z[2],X[3]>) ) / (Z[3] - X[3]):DDF := map(normal, DDF);Note that the (3,1) entry, in particular, is much simpler than the expression obtained in Example 1a.Verify that the defining equation for the NiMlJERERkc= matrix is satisfied.map(normal, f(X)-f(Z)-DDF.(X-Z));g1 := 'g1': g2 := 'g2': g3 := 'g3':<Text-field layout="Heading 2" style="Heading 2">Example 2</Text-field>An even starker contrast in the complexity of the formulas can be seen in the following example where NiMtJSJmRzYjJSJ5Rw== is a rational function.y := <y1, y2, y3>:
f := v -> <1, v[3], cos(v[1])/v[3]+(1-cos(v[1]))/v[2]>:
'f(y)' = f(y);Forming the Jacobian matrix and then applying straight-line integration to the (3,1) component, as in Example 1a, yields the following formula involving the functions NiMlI1NpRw== and NiMlI0NpRw== (the sine and cosine integrals).J := linalg[jacobian](f(y), y):J31 := unapply(J[3,1], y1,y2,y3);int(J31(X1+theta*(Z1-X1), X2+theta*(Z2-X2), X3+theta*(Z3-X3)), theta=0..1);In contrast, applying a dog-leg path of integration as in Example 1b yields the following simple entries in the DDF matrix.Rather than define each entry of the NiMlJERERkc= matrix by its explicit formula, as we did above, it is more compact to define the entries using procedure NiMlK3NjYWxhcmRpZmZH whose definition appears in the next section.n := Dimension(X):Xnew := copy(X):for j from 1 to n do
col[j] := <seq(scalardiff(f,i,j,Z[j],Xnew), i=1..n)>;
Xnew[j] := Z[j]
od:DDF := map(normal, linalg[augment](seq(col[j], j=1..n))):DDF := convert(DDF, Matrix);Note that the (3,1) entry, in particular, is vastly simpler than the expression obtained via straight-line integration.Verify that the defining equation for the NiMlJERERkc= matrix is satisfied.map(normal, f(X)-f(Z)-DDF.(X-Z));<Text-field layout="Heading 1" style="Heading 1">Algorithm for DDF Matrix Based on Scalar Divided Differences</Text-field>In [Geddes99] it is shown that the desired vector divided differences can be computed by the following method which relies on the computation of particular scalar divided differences.<Text-field layout="Heading 2" style="Heading 2">Procedure scalardiff</Text-field>The scalar divided differences required in the computation of a divided difference matrix NiMtJSREREZHNiQlIlhHJSJaRw== arise in the following context. The function NiMlImZH is a vector-valued function of a vector argument. It will be required to compute a scalar divided difference of NiMmJSJmRzYjJSJpRw== , a particular component of NiMlImZH , where we note that NiMmJSJmRzYjJSJpRw== is a scalar-valued function of a vector argument. Specifically, the following Maple procedure computes a scalar divided difference of NiMmJSJmRzYjJSJpRw== in the NiMlImpH-NiMlI3RoRw== variable, taking into account the case where it reduces to differentiation.scalardiff := proc(f::procedure, i::integer, j::integer, z::algebraic, X::Vector)
local g,upper,n,W,k;
# f maps an n-vector to an n-vector, where n = dim(X).
# If z <> X[j], the value returned is the scalar divided difference of
# g = f[i], the i-th component of f, in the j-th variable:
#
# ( g(X[1],...,X[j-1],z,X[j+1],...,X[n]) - g(X) ) / (z - X[j]) .
#
# If z = X[j], the value returned is the corresponding derivative value:
# g_j(X) where g_j = partial derivative of g w.r.t the j-th variable.
g := v -> f(v)[i];
if z <> X[j] then
upper := copy(X); upper[j] := z;
(g(upper) - g(X)) / (z - X[j])
else
n := Dimension(X); W := Vector(n, j -> w[j]);
eval(diff(g(W),W[j]), [seq(W[k]=X[k], k=1..n)])
fi
end:<Text-field layout="Heading 2" style="Heading 2">Theorem 1</Text-field>Let NiMlImZH be a vector-valued function of a vector argument where the vectors have dimension NiMlIm5H. Given any two NiMlIm5H-vectors NiMlIlhH and NiMlIlpH, define the divided difference matrix NiMtJSREREZHNiQlIlhHJSJaRw== byNiUvJi0lLH5+fn5+fn5+RERGRzYkJSJYRyUiWkc2JCUiaUclImpHLSUrc2NhbGFyZGlmZkc2JyUiZkdGK0YsJkYpNiNGLDcoJkYpNiMiIiIlJC4uLkcmRik2IywmRixGNkY2ISIiJkYoRjJGNyZGKDYjJSJuRy8lJn5+fn5pRztGNkY/LyUkfn5qR0ZCwhere NiMlK3NjYWxhcmRpZmZH is as defined above. Then the following identity holds:NiMvLCYtJSp+fn5+fn5+fmZHNiMlIlhHIiIiLSUiZkc2IyUiWkchIiIqJi0lJERERkc2JEYoRi1GKSwmRihGKUYtRi5GKQ== .<Text-field layout="Heading 2" style="Heading 2">Procedure DDFgen</Text-field>We can now present procedure NiMlJ0RERmdlbkc= which generates the matrix NiMtJSREREZHNiQlIlhHJSJaRw== defined in Theorem 1.DDFgen := proc(f::procedure, X::Vector, Z::Vector)
local i,j,n,Xnew,col;
# f maps an n-vector to an n-vector, where n = dim(X) = dim(Z).
# The value returned is an n-by-n matrix DDF(X,Z), computed via
# scalar divided differences, which satisfies the identity
# f(X)-f(Z) = DDF(X,Z) (X-Z) .
n := Dimension(X);
if Dimension(Z) <> n then ERROR("expecting dim(X)=dim(Z)") fi;
Xnew := copy(X);
for j from 1 to n do
col[j] := <seq(scalardiff(f,i,j,Z[j],Xnew), i=1..n)>;
Xnew[j] := Z[j]
od;
convert(map(normal, linalg[augment](seq(col[j], j=1..n))), Matrix)
end:<Text-field layout="Heading 1" style="Heading 1">A DDF-based Update Formula</Text-field>Recall the second order update formula based on the Jacobian:NiMvLSUiUUc2JCUiaEclImdHLChGKCIiIiomRidGKi0lImZHNiNGKEYqRioqKC0lIUc2IyomRiciIiNGNCEiIkYqLSYlIkpHNiNGLUYuRipGLEYqRio= .Consider two "nearby" vectors NiMlImdH and NiMlIlFH (which will correspond to iterates NiMmJSJ5RzYjJSJuRw== and NiMmJSJ5RzYjLCYlIm5HIiIiRihGKA==, respectively) and approximate the Jacobian NiMtJiUiSkc2IyUiZkc2IyUiZ0c= by NiMtJSREREZHNiQlImdHJSJRRw== . This leads to the following implicit update formula: Solve NiMvJSJRRywoJSJnRyIiIiomJSJoR0YnLSUiZkc2I0YmRidGJyooLSUhRzYjKiZGKSIiI0YyISIiRictJSREREZHNiRGJkYkRidGKkYnRic= for NiMvJSJRRy1GJDYkJSJoRyUiZ0c= .RemarksThis is an implicit formula, thus achieving desirable numerical stability properties.The computer algebra system makes it practical to form the DDF matrix for the given problem.In some cases of interest, it is possible to solve the above implicit formula for NiMlIlFH yielding an explicit symbolic solution!The latter point says that we can simultaneously achieve two apparently contradictory goals:the numerical stability properties of an implicit formula;the computational efficiency, in the purely numerical phase, achieved by an explicit formula for the next iterate, without the need to numerically solve a nonlinear system of equations at each step.<Text-field layout="Heading 1" style="Heading 1">Example 3</Text-field>This example illustrates the derivation of an explicit numerical update formula by solving implicit equations as discussed above.Consider the IVP specified by the following nonlinear second-order differential equation NiMlI2RlRw== with corresponding initial conditions specified by NiMlI2ljRw==.de := diff(x(t),t$2) - x(t)*diff(x(t),t) = 0;t0 := sqrt(2)*Pi/8: xt0 := 2-sqrt(2):ic := {x(t0) = xt0, D(x)(t0) = 2*xt0};For this example, the solution can be expressed in closed form: NiMvLSUieEc2IyUidEcqJi0lJXNxcnRHNiMiIiMiIiItJSR0YW5HNiMqJkYnRi1GKSEiIkYt .Converting NiMlI2RlRw== to a system of first-order equations, the initial-value problem becomesNiNDJC8tJSVkaWZmRzYkLSUieUc2IyUidEdGKy0lImZHNiNGKS8tJSh+fn5+fn55RzYjJkYrNiMiIiEmRilGNA==wherey := <y1,y2>:
f := v -> <v[2], v[1]*v[2]>:
'f(y)' = f(y);and the initial conditions are[y1,y2](t0) = [xt0, 2*xt0];Let us generate a numerical update formula by the NiMlJERERkc= approach of the previous section.g := <g1,g2>: Q := <Q1,Q2>:DDF := DDFgen(f,g,Q);Check the DDF defining property.map(normal, f(g)-f(Q)-DDF.(g-Q));The second-order implicit update formula of the preceding section is of the formNiMvJSJRRywoJSJnRyIiIiomJSJoR0YnLSUiZkc2I0YmRidGJyomLSUhRzYjKiZGKSIiI0YyISIiRictRi82IyooJSREREZHRiclIi5HRidGKkYnRidGJw==h := 'h':RHS := g + ScalarMultiply(f(g),h) + ScalarMultiply(DDF.f(g), h^2/2);Solving the implicit equations for the components of NiMlIlFH yieldssoln := solve({seq(Q[k]=RHS[k], k=1..2)}, {Q[1],Q[2]});solnlist := [subs(soln,Q1), subs(soln,Q2)];Q := evalf @ unapply(solnlist,h,g1,g2):Computing solution points via the explicit update formula is straightforward.h := 0.005;t0 := evalf(t0);g := [evalf(xt0), evalf(2*xt0)];L := [[t0,g[1]]]:for i from 1 to 330 do g := Q(h,g[1],g[2]); L := [op(L),[t0+i*h,g[1]]]od:plots[pointplot](L);The solution has a pole at NiMqJiUjUGlHIiIiLSUlc3FydEc2IyIiIyEiIg== = 2.221441469 . In this case, our numerical update formula has the traditional problem that it cannot integrate beyond a pole.<Text-field layout="Heading 1" style="Heading 1">Example 4: A Riccati Equation</Text-field>The following example illustrates that solving implicit equations for an explicit update formula can yield very interesting numerical methods. This example is discussed by Kahan [Kahan93] as an "unconventional numerical method". In this case, the DDF matrix is not used but rather a splitting technique is applied to the first-order update formula and then the implicit equations are solved.The IVP: NiNDJC8lI3knRywmJSJ0RyIiIiokJSJ5RyIiI0YoLy0lJn5+fn55RzYjIiIhRjA= .As an AIVP: NiNDJC8qJjckJiUieUc2IyIiIiZGKDYjIiIjRiolIidHRio3JEYqLCZGJ0YqKiRGK0YtRiovKiYlJX5+fn5HRiotRiY2JCIiIUY3Rio3JEY3Rjc= .g := 'g': h := 'h':y := <y1,y2>:
f := v -> <1, v[1]+v[2]^2>:
'f(y)' = f(y);g := <g1,g2>: Q := <Q1,Q2>:Recall that the first-order update formula (Euler's method) is NiMvLSUiUUc2JCUiaEclImdHLCZGKCIiIiomRidGKi0lImZHNiNGKEYqRio= . Apply a splitting technique to the second component of this vector equation (which corresponds to the original scalar IVP). The main idea is to replace the term NiMqJCUjZzJHIiIj by an average of NiMlI2cyRw== and NiMlI1EyRw==. Specifically, replace NiMtJSJmRzYjJSJnRw== in the first-order formula by NiMtJShmX3NwbGl0RzYlJSJoRyUiZ0clIlFH wheref_split := (h,v,w) -> <1, v[1] + h/2 + v[2]*w[2]>:f_split(h,g,Q);For curiousity, note that if one forms the NiMlJERERkc= matrix for this example one getsDDF := DDFgen(f,g,Q);but in what follows we are not using the NiMlJERERkc= matrix.Form and solve the implicit equations obtained via the splitting technique.RHS := g + ScalarMultiply(f_split(h,g,Q), h);soln := solve({seq(Q[k]=RHS[k], k=1..2)}, {Q[1],Q[2]}):solnlist := [subs(soln,Q1), subs(soln,Q2)]:Q := evalf @ unapply(solnlist,h,g1,g2);Applying this update formula yields the following results.h := 0.01;g := [0,0];L := [[g[1],g[2]]]:for step from 1 to 500 do g := Q(h,g[1],g[2]); L := [op(L),[g[1],g[2]]]od:plots[pointplot](L,view=[0..5,-50..50]);The true solution can be expressed in terms of quotients of Bessel functions. As pointed out by Kahan, one can verify that the above numerical method is performing admirably well through the poles of the solution, a property not common to conventional numerical methods.<Text-field layout="Heading 1" style="Heading 1">References</Text-field>[Geddes86] K.O. Geddes, Numerical integration in a symbolic context. Proceedings of SYMSAC'86, B.W. Char (ed.), ACM Press, New York, 1986, pp. 185-191.[Geddes92] K.O. Geddes and G.J. Fee, Hybrid symbolic-numeric integration in Maple. Proceedings of ISSAC'92, P.S. Wang (ed.), ACM Press, New York, 1992, pp. 36-41.[Geddes99] K.O. Geddes, Generating numerical ODE formulas via a symbolic calculus of divided differences. SIGSAM Bulletin 33 (128), 1999, pp. 29-42.[GeddesGonnet89] K.O. Geddes and G.H. Gonnet, A new algorithm for computing symbolic limits using hierarchical series. Appears in Symbolic and Algebraic Computation, P. Gianni (ed.), Lecture Notes in Computer Science, No. 358, Springer-Verlag, Berlin, 1989, pp. 490-495.[Kahan93] W. Kahan, Unconventional numerical methods for trajectory calculations. Manuscript, prepared for the Workshop on Integration Algorithms for Classical Mechanics, Fields Institute for Research in Mathematical Sciences, Waterloo, Ontario, Canada, Oct. 1993.[Kahan97a] W. Kahan and Ren-Cang Li, Unconventional schemes for a class of ordinary differential equations -- with applications to the Korteweg-de Vries equation. J. Comp. Physics 134, 1997, pp. 316-331.[Kahan97b] W. Kahan and Ren-Cang Li, Composition constants for raising the orders of unconventional schemes for ordinary differential equations. Math. Comp. 66, 1997, pp. 1089-1099.[Richardson96] D. Richardson, B. Salvy, J. Shackell and J. Van derHoeven, Asymptotic expansions of exp-log functions. Proceedings of ISSAC'96, Y.N. Lakshman (ed.), ACM Press, New York, 1996, pp. 309-313.[Salvy91] B. Salvy, Examples of automatic asymptotic expansions. SIGSAM Bulletin 25 (2), 1991, pp. 4-17.[Salvy92] B. Salvy, General asymptotic scales and computer algebra. Appears in Asymptotic and Numerical Methods for Partial Differential Equations, Critical Parameters and Domain Decomposition, H. Kaper and M. Garbey (ed.), Kluwer Academic Publishers, 1992, pp. 255-266.TTdSMApJM1JUQUJMRV9TQVZFLzQzMDcyOFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCIiIiUjeTNHLCYqJi0lJGNvc0c2IyUjeTFHRidGCihGJ0YnKiYsJkYnRidGKyEiIkYnJSN5MkdGJ0YnNiIKTTdSMApJM1JUQUJMRV9TQVZFLzUyMTUwMFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIhRicqJiwqKiYtJSRjb3NHNiMlI1oxRyIiCiIlI1gzR0YvISIiKiYlI1gyR0YvRitGL0YvKiYtRiw2IyUjWDFHRi9GMEYvRi8qJkYzRi9GNUYvRjFGLywmRi5GMUY3Ri9GMUYnRicsJkYKL0YvRitGMUYnRi9GKzYiCg==TTdSMApJM1JUQUJMRV9TQVZFLzQzMDgwOFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCIiIUYnRic2Igo=TTdSMApJM1JUQUJMRV9TQVZFLzUzMDMyMFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCIiIiUjeTNHLCYqJi0lJGNvc0c2IyUjeTFHRidGCighIiJGJyomLCZGJ0YnRitGL0YnJSN5MkdGL0YnNiIKTTdSMApJM1JUQUJMRV9TQVZFLzUzNzIwOFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyoiJCIkIiIhRicqKiwqKiYlI1gyRyIiIi0lJGNvc0c2CiMlI1oxR0YsISIiKiZGLUYsJSNYM0dGLEYsKiZGK0YsLUYuNiMlI1gxR0YsRiwqJkY1RixGM0YsRjFGLEYrRjFGM0YxLCZGMEYxRjdGLEYKMUYnRicqKCwmRjFGLEYtRixGLCUjWjJHRjFGK0YxRidGLCwkKihGLUYsJSNaM0dGMUYzRjFGMTYiCg==TTdSMApJM1JUQUJMRV9TQVZFLzUzMDQ0MFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiQiJCIiIUYnRic2Igo=TTdSMApJM1JUQUJMRV9TQVZFLzUzMDQ4MFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiMiIyUjeTJHKiYlI3kxRyIiIkYnRio2Igo=TTdSMApJM1JUQUJMRV9TQVZFLzU1MTU2MFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyUiIyIjIiIhJSNnMkciIiIlI1ExRzYiCg==TTdSMApJM1JUQUJMRV9TQVZFLzUzMDUyMFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiMiIyIiIUYnNiIKTTdSMApJM1JUQUJMRV9TQVZFLzUzMDU2MFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiMiIywoJSNnMUciIiIqJiUiaEdGKSUjZzJHRilGKSooCkYrIiIjRihGKUYsRikjRilGLiwoRixGKSooRitGKUYoRilGLEYpRikqJkYrRi4sJiokRixGLkYpKiglI1ExR0YpRihGKUYsRilGKUYpRi8KNiIKTTdSMApJM1JUQUJMRV9TQVZFLzUzMDYwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiMiIyIiIiwmJSN5MUdGJyokJSN5MkciIiNGJzYiCg==TTdSMApJM1JUQUJMRV9TQVZFLzUzMDY0MFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiMiIyIiIiwoJSNnMUdGJyUiaEcjRiciIiMqJiUjZzJHCkYnJSNRMkdGJ0YnNiIKTTdSMApJM1JUQUJMRV9TQVZFLzU3NjE0OFgsJSlhbnl0aGluZ0c2IjYiW2dsISIlISEhIyUiIyIjIiIhIiIiRicsJiUjUTJHRiglI2cyR0YoNiIKTTdSMApJM1JUQUJMRV9TQVZFLzUzMDY4MFgqJSlhbnl0aGluZ0c2IjYiW2dsISMlISEhIiMiIywmJSNnMUciIiIlImhHRiksJiUjZzJHRikqJkYqCkYpLChGKEYpRiojRikiIiMqJkYsRiklI1EyR0YpRilGKUYpNiIK