Some Russian OCR'ed the article on the web. It's pretty trashed, but you might get something out of the text.
http://booklists.narod.ru/_Papers/Computer_algebra/Stoutemeyer._Solving_integral_equations_in_CAS_T__19s_.1.htm
At 02:45 PM 7/14/2011, Raymond Toy wrote:
>>>>>> "John" == John Ogilvie <ogilvie at cecm.sfu.ca> writes:
> John> Does anybody have a copy of a report by David R. Stoutemyer of title
> John> "Analytical solution of integral equations using computer algebra"
> John> issued as Computer Physics Reports No. 34, University of Utah, Salt Lake
> John> City, Utah USA, 1975?
> John> I should be grateful for a copy of this report.
> John> John Ogilvie
>
>Isn't this usually solved by visiting your local university library?
>Or, at worst, interlibrary loan?
>
>Sorry, I don't have this report.
>
>Ray
Stoutemeyer. Solving integral equations in CAS(T)(19s)
Top / Papers / Computer algebra / [1]
--page0001--
Analytically Solving Integral Equations
by Using Computer Algebra
DAVID R. STOUTEMYER
University of Hawaii
This paper describes how a computer algebra language such as REDUCE can be used to
automatically construct closed-form and series analytical solutions of integral equations.
The integral equations may be from a large class of nonsingular scalar one-dimensional
problems, including nonlinear generalizations of the classical Volterra and Fredholm equations.
Given a symbolic expression for the integral equation, the program described here categorizes
the equation and applies applicable symbolic techniques to construct the solution (s). If no
closed-form solutions are achieved, then series solutions are automatically attempted; and
if they too are unsuccessful, then approximate solutions of the collocation type are attempted.
In any case, the solutions differ from numerical floating-point solutions in the following
respects.
A) The solution expression is a complete direct representation of the solution for all values
of the independent variable over an interval rather than a direct representation of the solution
at a finite number of points or an indirect representation consisting of the numerical coefficients
intended for use in an expression or interpolation formula.
B) The equation may have nonnumerical parameters such as symbolic integration limits,
scale factors, or coefficients, with the solution being expressed in terms of these parameters.
C) Any arithmetic is done exactly by using arbitrary-precision rational numbers
Key Words and Phrases: integral equations, Volterra equations, Fredholm equations,
REDUCE, symbolic algebraic manipulation, polyalgorithm
CR Categories: 5.7, 5.15, 5.18
1. INTRODUCTION
Throughout applied mathematics, closed-form analytical solutions are particularly
desirable because of the insight they afford. Such solutions can provide a complete
representation over the entire range of data, and, provided they are not too
complicated, such solutions permit direct qualitative understanding of how the solution
varies with the data. Unfortunately this kind of solution is rarely possible except
for simple cases that arise only in classrooms or from drastic idealizations of real
physical situations. Consequently we must often be content with a less informative
numerical solution, perhaps together with a numerical sensitivity study. For many
Copyright ? 1977, Association for Computing Machinery, Inc. General permission to
publish, but not for profit, all or part of this material is granted provided that ACM's
copyright notice is given and that reference is made to the publication, to its date of issue, and
to the fact that reprinting privileges were granted by permission of the Association for
Computing Machinery.
This work was supported in part by the National Science Foundation under Grants GJ-32181
and MCS75-22983 and in part by the Advanced Research Projects Agency, Office of the
Department of Defense under Contract DAHC 15-73-C-0363.
Authors' address: Department of Electrical Engineering, University of Hawaii, Honolulu,
Hawaii 96822.
ACM Transactions on Mathematical Software, Vol. 3, No 2, June 1977, Pages 128-146
.
--page0002--
Analytically Solving Integral Equations ? 129
areas of applied mathematics, however, computer-algebra languages permit
extension of the threshold of difficulty beyond which analytical techniques must be
abandoned.
The solution of various kinds of equations is one of the major endeavors of
applied mathematics, and general-purpose programs have already been written for
solving equations in finite terms and differential equations. Examples are the
programs of Yun [23], Moses [14, 15], Golden and Kuipers [7], Bogen [2], and Kulp
[10]. Integral equations are probably the next most frequently encountered
category in engineering and science, and computer-algebra programs have been written
to solve a particular integral equation (Gentlemen [6]) and to implement one
technique for solving a certain class of integral equations (the Mathlab group
[11])-
This paper describes a general-purpose program intended to construct
form series, or projected analytical solutions to a large class of one-dimensional
integral equations. The program utilizes many of the most widely applicable
methods for constructing such solutions. Section 2 outlines these methods, and Section 3
presents elementary examples using the corresponding program. Section 4
summarizes test results for a variety of more challenging problems, with conclusions
given in Section 5.
2. IMPLEMENTED ANALYTICAL TECHNIQUES
The type of integral equation considered here is of the form
Q B)
I. Ja(x) )
and/or of the first kind
I w[x,u,p(u)]du = f(x) C)
Ja{x)
that are collectively equivalent to eq. A). For example, the equation
\" If fb \\ ' 1
in / p(u) du log< p(x) ? / exp [p(u)] du> ? 2 asin I p(u) du
asm
- / exp [p(u)]du> + 2 =
ACM Transactions on Mathematical Software, Vol. 3, No 2, June 1977
.
--page0003--
130 ? D. R. Stoutemyer
would be converted to the pair of independent integral equations
I p(u) du = sin (?;), D)
J?x
p(x) = / exp [p(u)] du + e. E)
SOLVE, described by Stoutemyer [18], recursively uses factorization, together
with known inverses of various polynomial, algebraic, and transcendental functions,
to solve a nonlinear equation in finite terms. SOLVE is also capable of symbolically
solving simultaneous algebraic equations provided they are linear or linear
fractional in the unknowns.
Most of the solution techniques require a specialized form of eqs. B) or C),
such as the second kind
p(x)=y F)
I J?fa) J
or quasi-linear second kind
p(x) = J w[x,u,p(u)]du+f(x) G)
or a linear equation, for which w in eq. C) or G) is the form k(x, u)p{u). Also,
some techniques are applicable only for fixed limits, where the limits of integration
are free of x, and others are applicable only for variable limits, where one or both
limits of integration depend on x. Usually at most one limit is variable, and it is
simply a; or ? x; but it is easier to program the general case than to transform the
usual cases to one standard form and check for violations.
The integral-equation solver automatically checks the necessary consistency
condition for variable-limit equations of the first kind
/(?) =0, Yx b(f) = a(x). (8)
For each of the (consistent) independent simplified equations, the program
successively attempts appropriate techniques until one succeeds or none remain.
Depending upon the particular form of the simplified integral equation, the techniques
are attempted in approximately the order that they are discussed below. Partial
successes, such as a closed-form solution containing an unevaluated integral or a
series solution with less than the requested number of terms, are retained too, but
without terminating the search. The user also has the option of specifying that
all appropriate techniques be attempted regardless of any successes, or specifying
that only a particular technique be attempted. Originally intended as debugging
tools, these options increase the flexibility of the program by facilitating survey
and comparison of the techniques, followed perhaps by a concentrated effort on
the most promising technique(s)?one(s) which otherwise might have been
preempted by a previous success or by exhaustion of the available storage while a
previous technique was being attempted.
Table I summarizes the individual techniques described herein.
ACM Transactions on Mathematical Software, Vol 3, No 2, June 1OT.
.
--page0004--
Analytically Solving Integral Equations * 131
Table I. Summary of Implemented Techniques
Method
Fixed-limit finite-rank 2nd kind
Variable-limit rank-1 2nd kind
Variable-limit rank-1 1st kind
Laplace transform
Fourier transform
Differentiation
Fredholm-Carleman series
Taylor series
Neumann series
1st kind series
Collocation
Selector
FINITERANK
FINITERANK
FINITERANK
TRANSFORM
TRANSFORM
?
FREDSERIES
TAYLOR
NEUMANN
FIRSTICINDSERIES
COLLOCATE
Equation factor convertible
to the form:
p(x) = y{x,\\ I q.(x)r {u,p(u)\"\\du}
J jml 3 0
b(x)*
p(x) = f(x) + j q(x)R(u)p(u)du
aix)
bix)*
f q(x)r[u,p(u)]du = fix)
a(x)
X
[0 or p(x)] = f(x) + \\K(x-u)p(u)du
0
oo
[0 or p(x)J = f(x) + \\K(x-u)p(u)du
?-co
b(x)*
[ [ , , = fix)
aix)
bix)
pix) =f(x) + K(x,u)p(u)du
aix)
bix)*
pCx) =fix) + j U{x,u,p(u)]du
aix)
pix) = { , [ ( )] ,
bix)
| U[x,u,pt$iu,x))*ldu}
aix)
bix)
\\w{x,u,p[$(u,x)]f}du = fix)
aix)
either of the above two forms
* Either aix) ' = 0 and bix) =x or bix) ' = 0 and aix) = x,
f There may be any number of such arguments, each with a distinct a or B,
including the usual degenerate cases a(x)=x, Biu,x)=u.
ACM Transactions on Mathematical Software, Vol. 3, No. 2, June 1977.
.
--page0005--
132 ? D. R. Stoutemyer
2.1. Fixed-Limit Finite-Rank Equations
An integrand in F) of the form
?^ , (?I (9)
3-1
is called finite-rank, degenerate, or separable.
For fixed limits substitute eq. (9) into eq. F), then interchange the order of
summation and integration to give
p(x) = .}^)], A0)
where the unknown constants c, are given by
fb
c, = / r,[u, p(u)] du, j = 1, 2, . . ., n. A1)
\"a
Substituting for x and i for j in eq. A0) gives
!>.?.(?)]. A2)
Substituting eq. A2) into eq. A1) gives
c3 = I rAu, \\u, ^ . ( ) \\>du, j=l,2,...,n, A3)
which is a set of n simultaneous equations in finite terms for C\\ through cn . After
these equations have been solved by using symbolic integration and the SOLVE
function, substitution of ci through cn into eq. A0) gives the exact closed-form
solution. As mentioned before, SOLVE is currently limited to linear or
fractional equations for n > 1, but there are plans to incorporate Yun's algorithm
for simultaneous polynomial equations [23]. (Davis [4] and Squire [21, pp.
233], solve some nonlinear finite-rank examples by hand.) Meanwhile, the
equation solver prints eq. A0) together with eq. A3) when SOLVE cannot solve
eq. A3), permitting the user to proceed numerically with a problem that is usually
easier than a direct numerical treatment of the original integral equation. The
integral-equation solver is written in the REDUCE computer-algebra language,
which does not currently have a built-in integrator. Consequently the symbolic
integrations are performed by a simple REDUCE program described by
Stoutemyer [19].
By using factorization and collection of terms together with substitutions such
aslog[a(a:)|fl(u)]?>loga(x) + log/?(w), sin [a(x) + @(u)] ?> sin [a(x)] cos \\0(u)]
+ sin \\j3(u)} cos [a(x)}, and c[aM+Mu)} -? calx)c0w, the program attempts to
minimize n and to achieve the form of eq. (9) even when w is not originally expressed
that way.
2.2 Rank-1 Second-Kind Variable-Limit Equations
As developed by Cochran [3], the general closed-form solution to the linear rank-1
variable-limit equation of the second kind, with w = q(x)R(u)p(u), is
ACM Transactions on Mathematical Software, Vol 3, No 2, June 1977
.
--page0006--
Analytically Solving Integral Equations ? 133
p(x) = f(x) + q(x) [ U> f{u)R{u) exp \\f q(v)R(v) dv~\\du, A4)
provided either b(x) ? x and d{x) = 0 or a(x) = a; and b'{x) = 0, as is usual.
2.3 Rank-1 First-Kind Variable-Limit Equations
For a rank-1 variable-limit equation of the first kind, assume sufficient
differentiability and q(x) ?? 0 for a(x) < b(x), with w = q(x)r[u, p(u)]. Then division
followed by differentiation gives
{ ( ), p[b(x)]}b'(x) - rla(x), p[a(x)])a'(x) = U(x)/q(x)}', A5)
where primes denote differentiation with respect to x. This is an equation in finite
terms for p(x), provided either b(x) = a;anda'(a;) = ( ) = x&ndb (x) = 0,
as is usual.
2.4 Integral Transforms
Laplace and Fourier transforms provide closed-form general solutions to certain
linear equations with convolution or difference type kernels,
k(x, ) = { ? ), A6)
which may be checked by determining if k(x, + ) is free of x. We can use the
Laplace transform for variable-limit equations when the lower integration limit is
0 and the upper limit is x, and we can use the Fourier transform when the lower
limit is ? oo and the upper limit is ??. Let denote the appropriate transform of
and F denote the appropriate transform of/. Then the solution to the first kind
of equation is the inverse transform of
F/K, A7)
and the solution to the second kind of equation is the inverse transform of
F/(l - K). A8)
REDUCE does not have built-in Laplace or Fourier transforms; so a modest
transform program listed in Stoutemyer [20] has been written for testing purposes
to serve until completion of a more powerful program under current development.
2.5 Differentiation
Provided there is sufficient differentiability, a variable-limit equation of the first
kind can be differentiated to give
w{x, b(x), p[b(x)]}b'(x) - w{x, a(x), p[a(x)]]d(x)
wx[x,u,p(u)]du-f(x)=O, A9)
where the x subscript denotes partial differentiation with respect to x. This may be
an equation of the second kind, provided either b(x) = x and a (x) = 0 or a(x) = x
and b (x) = 0, as is usual. The integral-equation solver is then recursively reentered
with this new equation. If eq. A9) has factors of the first kind, such as when the
first two terms of A9) are identically zero, the process may be reapplied to these
ACM Transactions on Mathematical Software, Vol. 3, No. 2, June 1977
.
--page0007--
134 ? D. R. Stoutemyer
factors. However, the depth of recursion is limited to a user-specified value to
prevent infinite looping.
2.6 Fredholm-Carlemdn Series
For linear integral equations of the second kind, there are well-known recurrence
relations for constructing the Fredholm series solution p. Beginning with g = G
= k(x,u) andj = D = 1, successively iterate the following assignments:
Mx)
/?OI.X)
d <- / - g(u, u) du/j,
Ja(x)
/?biz)
g <? k(x, u)d + / k(x, v)g(v, a) dv,
D<^-d + D,
G<-g + G,
Mx)
f> *-fix) + / G(x, u)f(u) du/D.
Although the series terminates after a finite number of terms for the fixed-limit
finite-rank case, the program uses the method of Section 2.1 when it can do so.
2.7 Taylor Series
Noble [16] describes a Taylor-series technique that is applicable to quasi-linear
variable-limit equations of the second kind, and Norman [17] uses algebraic
manipulation to apply it to a particular integral equation. Let x denote a solution of
a(x) = b(x) and assume p(x) can be represented in the form
P(x) = /(x) + g {X T, *Y ~ b>{x) - f(x)]x=i B0)
for sufficiently small | x ? x \\. A violation of this assumption, such as a singularity
at x, in lp(x) ? f(x)] or one of its derivatives is generally revealed by a
divide interrupt during construction of the series. The program traps these errors
from this source and prints the more specific message: \"Zero-divide when
constructing Taylor series. Probable cause is a singularity in the solution or one of its
derivatives.\" From eq. G),
lp(x) - /(*)] = 0,
b(*) -f(x)]Li = {w[x, b(x), p(b(x))]b'(x) - w[x, a(x), p(a(x))]a'(*)}_s
= w[x,x,f(x)][b'(x) -a'(a)]x=i,
provided a(x) = b(x) = x, as is usual.
The Taylor-series method, when applicable, has the advantage over the
holm-Carleman and Neumann series of usiiig symbolic differentiation rather than
the more difficult process of symbolic integration.
2.8 NeumaHh Series
The Neumann series is applicable to equations of the quasi-second kind B). It is
in fact applicable to the more general situation where one or more subexpressions
ACM Transactions on Mathematical Software, Vol 3, No 2, June 1977
.
--page0008--
Analytically Solving Integral Equations ? 135
of the form p[a(x)] occur in 9, where one or more subexpressions of the form
p\\0(x, u)] appear in the integrand, and when there is more than one integral. This
method is simply successive approximation: It begins with a guess for p that is
substituted into the right side of eq. B), which is then evaluated to provide a new
approximation to p. This new approximation is then substituted into the
hand side, and so on. This process may of course diverge if the right-hand side of
eq. B) is too sensitive to its second or third argument. Often the integrand
includes a scale factor, and the series converges only for sufficiently small magnitudes
of that factor. The form of successive partial sums may indicate the critical scale
factor; but alone they do not constitute a proof that it is critical.
The user may specify the first guess, which may influence the success of the
symbolic integrations, but the default first guess is the right side of eq. B) with
the integrals set to zero.
Fortuitously exact results are recognized a.nd are so indicated by the program.
2.9 First-Kind Series
Mikhlin and Smolitskiy [12] give a series which is sometimes successful for equations
of the first kind: Beginning with an arbitrary guess p* and a coefficient X, iterate
using
x) - / w[x, u, f(u)]> du. B1)
?/--page0018--
Analytically Solving Integral Equations ? 145
boolean procedure RANKlSECOND(y, P_OF_X, IESOLN);
begin comment Equation is variable-limit of the form P_OF_X=y.
Determines if is linear. If so, uses SUMFACTORS to determine if the
integrand is convertible to the form q(x)R(u)p(u). If so, uses equation
04) and appends result together with \"FINITERANK\" on IESOLN. Returns
true or false according to the success of its endeavor;
end RANKlSECOND;
boolean procedure RANK1FIRST(F, P_OF_X, INTGRL, IESOLN);
begin comment Equation is variable-limit of the form INTGRL=F,
with INTGRL of the form INTEGRAL(w, u, a(x), b(x)). Uses SUMFACTORS to
determine if w has the form of equation (9) with n=l. If so, SOLVE is
used on equation A5), and solutions paired with \"FINITERANK\" are appended
to IESOLN. Returns true or false according to the success of its endeavor;
... end RANK!FIRST;
symbolic procedure SUMFACTORS (w, x, u);
begin comment If w can be put in the.form of equation (9), returns a list
of the form [(qi>n), D2>r2)? ???? ( > )]- Returns empty otherwise;
symbolic q0, r0;
w -<- FR(w, x, u);
if w contains any FR(...) then return empty;
Factor the denominator of w, collecting into q0 the product of all factors
that are free of u, and collecting into r0 the product of all remaining
factors that are free of x, If any factors remain then return empty;
4o \"~ Vqo'. r0 +? l/r0;
Factor the numerator of w, collecting into qo the product of all factors that
are free of u, and collecting into r0 the product of all remaining factors
that are free of x;
Expand the product of any remaining factors of w with the x and all functional
forms of x ordered before u and all functional forms of u to get a form
such as (9) with independent monomial qj expressions. Repeat for the
reverse ordering to get a form such as (9) with independent monomial rj
expressions, and use the second expansion if it gives a smaller value of n
than the first expansion;
Multiply each of the qj by q0 and multiply each of the rj by r0, then return
a list of these qj r,- pairs;
end SUMFACTORS; J
comment pattern-matching rules for the operator FR, intended to convert some
transcendental expressions into finite-rank;
for all Y, Z, X, U let FR(-Y,X,U) + -FR(Y,X,U),
FR(Y-Z,X,U) + FR(Y,X,U)-FR(Z,X,U),
FR(Y+Z,X,U) + FR(Y,X,U)+FR(Z,X,U),
FR(Y*Z,X,U) ->- FR(Y,X,U)*FR(Z,X,U),
FR(Y/Z,X,U) + FR(Y,X,U)/FR(Z,X,U);
if free of (arbitrary X, arbitrary R_OR_Q) or free of (arbitrary U, R_OR_Q) let
FR(R_OR_Q?X,U) -, R_or_Q;
if free of (arbitrary X, arbitrary R) and free of (arbitrary U, arbitrary Q)
and not free of (U,R) and not free of (X,Q) let
FR(SIN(.Q+R),X,U) ->- SIN(Q)*COS(R)+COS(Q)*SIN(R),
FR(COS(Q+R),X,U) + COS(Q)*COS(R)-SIN(Q)*SIN(R),
FR(LOG(Q/R),X,U) ?* LOG(Q)-LOG(R),
FR(LOG(R/Q),X,U) ->- LOG(R)-LOG(Q),
FRCLOG(Q*R),X,U) + LOGCQ)+LOG(R);
boolean procedure SOLVEANDSUB (EXPRESSIONS, UNKNOWNS, FORM, TECHNIQUE, IESOLN);
begin comment EXPRESSIONS is a list of one or more expressions and
UNKNOWNS is a list of unknowns. TECHNIQUE is either \"FINITERANK\" or
\"COLLOCATE\". FORM is an expression containing the indeterminates listed
ACM Transactions on Mathematical Software, Vol. 3, No. 2, June 1977.
.
--page0019--
146 ? D. R. Stoutemyer
in UNKNOWNS. The equation(s) EXPRESSIONS=O are solved for UNKNOWNS.
For each solution, the values of UNKNOWNS are substituted in FORM, and the
result paired with TECHNIQUE is appended to IESOLN. Returns true or false
according to the success of this endeavor;
... end SOLVEANDSUB;
boolean procedure TRANSFORM (y_0R_F, P_QF_X, INTGRL, IESOLN).
begin comment If INTGRL=O then equation is of the 2nd kind, with P_OF_X=y_OR_F.
Otherwise it is of the 1st kind with INTGRL=y_OR_F, If the equation is of
linear convolution type with appropriate integration limits, uses formula A7)
or A8) with Laplace or Fourier transforms depending upon the kind and the
integration limits. Appends any solution together with \"TRANSFORM\" onto
IESOLN and returns true or false according to the success of this endeavor.
... end TRANSFORM;
boolean procedure FREDSERIES (y, P_OF X, N, IESOLN);
begin comment Equation is fixed-limit, of the form P_OF_X=y. N is the
requested number of terms. If is nonlinear, then false is returned
immediately. Otherwise, the iteration of section 2,6 is executed up to
N times, or until one of the Integrals cannot be evaluated in closed form,
p from the last iteration, together with \"FREDSERIES\" is appended to IESOLN.
Returns true or false according to the success in achieving the requested
number of iterations without unevaluated integrals;
... end FREDSERIES;
boolean procedure TAYLOR Cy, P_OF X, N, IESOLN);
begin comment Equation is variable-limit of the form P_OF_X=y. N is the
requested number of terms. Determines if is quasi-linear. If not then
false is returned immediately. Otherwise, the iteration of section 2.7 is
repeated N times, or until a singularity in the solution or one of its
derivatives causes a zero-divide interrupt. Appends the series together
with \"TAYLOR\" onto IESOLN and returns true or false according to the
success of the endeavor;
... end TAYLOR;
boolean procedure FIRST_OR_NEUMANN (y_0R_F, P_0F_X, INTGRL, N, GUESS, IESOLN);
begin comment If INTGRL^O then the equation is of the form rNTGRL=y_OR_F.
Otherwise the equation is of the form P_OF_X=y OR_F. Accordingly either
the First-kind iteration of section 2.9 or the~~Neumann iteration of section
2.8 is repeated up to times, beginning with GUESS, until an integral
cannot be evaluated in closed form. Appends any solution together with
appropriately \"FIRSTKINDSERIES\" or \"NEUMANN\" to IESOLN, then returns
true or false according to the success of the endeavor;
... end FIRST_OR NEUMANN;
boolean procedure COLLOCATE (EXPRESSION, P_OF_X, N, IESOLN);
begin comment Equation is of the form EXPRESSIONS; P_OF X is of the form P(x).
Substitutes APPROX(..,,N) for every instance of P(,.J in EXPRESSION,
then evaluates EXPRESSIONS at POINT(a,b,l ,N), POINT,N),...,
POINT(a,b,N,N) to get a set of N equations in terms of the N indeterminates
in the expression returned by APPROX, If none of the N equations contains
an unevaluated integral, then SOLVEANDSUB is used to solve the equations
for the N indeterminates and substitute them into APPROX(X.N). Prints the
corresponding equations and APPROX(X.N) if SOLVEANDSUB fails. Otherwise
SOLVEANDSUB appends the solutions to IESOLN. Returns true or false
according to the success of the endeavor; ... end COLLATE;
symbolic procedure APPROX(X, N);
comment The user may wish to replace this default definition of APPROX;
symbolic procedure POINT(A, B, J, N);
comment The user may wish to replace this default definition of POINT;
if N=1 then (A+B)/2
else A+(B-A)*(J-l)/(n-l);
ACM Transactions on Mathematical Software, Vol 3, No. 2, June 1977.