Compact Finite Difference Method for Integro ... - Semantic Scholar
Compact Finite Difference Method for IntegroDifferential Equations Jichao Zhao and Robert M. Corless Abstract. In this paper, we give sixth order compact finite difference formula for second order integro-differential equations (IDE) with different boundary conditions, and both of error estimates and numerical experiments confirm our compact finite difference method can get fifth order of accuracy. We also adjust compact finite difference method for first order IDE and a system of IDE and give numerical experiments for them. Our algorithm even can solve nonlinear IDE and unsplit kernel of IDE. The most advantages of compact finite difference method for IDE are that it obtains high order of accuracy, while the time complexity to solve the matrix equations after we use compact finite difference method on IDE is O(N ), and it can solve very general case of IDE. Keywords. Compact finite difference method, IDE, Integro-differential equations, Fredholm equations, Volterra equations, High accuracy.
1. Introduction Integro-differential equation (IDE) is an equation that the unknown function appears under the sign of integration and it also contains the derivatives of the unknown function. It can be classified into Fredholm equations and Volterra equations. The upper bound of the region for integral part of Volterra type is variable, while it is a fixed number for that of Fredholm type. For example, the following integro-differential equation belongs to Volterra type: Z x d2 u(x) + K(x, s)u(s)ds = f (x, u), a ≤ x ≤ b (1.1) − dx2 a and equation (1.2) is Fredholm integro-differential equation. In this paper we focus on Fredholm integro-differential equations, but all algorithms in our work can be applied to Volterra integro-differential equation with a little modification. If the kernel function holds K(x, s, u(s)) = K(x, s)u(s), we call the kernel is “split” type,
2
Jichao Zhao and Robert M. Corless
otherwise it is “unsplit” type. Our compact finite difference method works even for unsplit kernel, to simplify, we write the kernel function separately. Integro-differential equations arise from many branches of science, for example in control theory and financial mathematics (refer to [15], [8], [13], and [14]). Integro-differential equations are important, but they are hard to solve even numerically, so the progress on how to solve them is slow. In this paper, we mainly deal with the following integro-differential equation [6]: Z b d2 u(x) − + K(x, s)u(s)ds = f (x, u), a ≤ x ≤ b, (1.2) dx2 a where function f (x, u) can be a nonlinear function of u(x), and the typical boundary conditions: u(a) = ua , u(b) = ub ,
(1.3)
here ua and ub are already known boundary values. We assume the above equation (1.2) under the following two assumptions [6]: 1. The forcing function f (x, u) and all its partial derivatives are continuous , (x,u) and ∂f ∂u is nonpositive for a ≤ x ≤ b. 2. The kernel function K(x, s) satisfies the following positive definite property: Z bZ b K(x, s)ξ(x)ξ(s)dsdx > 0, (1.4) a
a
where ξ(x) is any continuous nonzero function, and holds Z bZ b | K(x, s) |2 dsdx < +∞. a
(1.5)
a
We know from the references [2], [6], and [14], the problem (1.2∼1.3) has a unique solution u(x) under the above two assumptions. In section 2, we give a brief introduction to a high accurate compact finite difference formula for second derivatives. In section 3, we show how to use compact finite difference method to solve the second order integro-differential equations with varying boundary conditions, also give the algorithm and error estimates for them. Numerical experiments in the last part of this section show that compact finite difference method can obtain fifth order of accuracy. In section 4, we give a high accurate compact finite difference method for first derivatives and use it to solve first order integro-differential equations, then followed by numerical examples to show this method also obtains fifth order of accuracy. As an expansion, we describe how to solve a system of integro-differential equations and an illustration is given in section 5. Conclusions are drawn in section 6. Note we use v(x) and u(x) to denote the numerical solution and exact solution of integro-differential equations respectively in the following sections.
Compact Finite Difference Method for IDE
3
2. Sixth Order Compact Finite Difference Formula Compact finite difference method is a special finite difference method which uses the values of the function and its derivatives only at three consecutive points. First discretize the interval [a, b], let spatial step h = 4x = |b−a| N , xi = a + ih, i = 0, 1, · · · , N − 1, N . 2.1. Standard Compact Finite Difference Formula The standard sixth order compact finite difference formula for second derivatives is d2 vi d2 vi+1 h2 d2 vi−1 ( + 10 + ) = vi−1 − 2vi + vi+1 , (2.1) 12 dx2 dx2 dx2 where vi = v(xi ), and the coefficients can be determined in the following way: 1. Write the desired compact finite difference formula with unknown coefficients: d2 vi−1 d2 vi d2 vi+1 + m0 2 + m1 ) = a−1 vi−1 + a0 vi + a1 vi+1 , (2.2) 2 dx dx dx2 here m−1 , m0 , m1 , a−1 , a0 , and a1 are the parameters to be decided.
h2 (m−1
2. Expand both sides of the equation (2.2) by Taylor series at the point xi with respect to the discretization parameter h, then collect them by the order of h. 3. We can obtain six equations by setting the coefficients of hj , j = 0, 1, · · · , 4, 5, equal zero. Finally solve the six equations for the six unknown parameters. From the above algorithm, we can see easily the accuracy is O(h6 ) for formula (2.1). For more details on how to generate compact finite difference formula, please refer to [7], [11], and [12]. Next we will show how to adjust compact finite difference formula for differential equations. 2.2. Solving Differential Equations by Compact Finite Difference Method Suppose we want to solve the following second order differential equation numerically for u(x): d2 u(x) = f (x), a ≤ x ≤ b, with u(a) = ua , u(b) = ub , (2.3) dx2 where f (x) is a known function. Write equation (2.3) in a discrete point form: −
d2 u(xi ) = f (xi ), 1 ≤ i ≤ N − 1, with u0 = ua , uN = ub . (2.4) dx2 Now the problem is that we know the values of second derivatives u00 (xi ) by equation (2.4) and need to solve for u(xi ), 1 ≤ i ≤ N − 1. So if we can find some relation between u00 (xi ) and u(xi ), then we can solve the equation (2.3) numerically. The tradition way is to use the second accurate central difference formula, while compact finite difference formula provides a higher accurate method to solve equation (2.3). We use standard sixth order compact finite difference formula (2.1), −
4
Jichao Zhao and Robert M. Corless
i = 2, · · · , (N − 2), for interior points. Since we have known the boundary conditions, we need to adjust compact finite difference formula for the boundary points. When i = 1, we use h2 (14v100 − 5v200 + 4v300 − v400 ) = v0 − 2v1 + v2 , 12 and when i = N − 1, h2 00 00 00 00 (−vN −4 + 4vN −3 − 5vN −2 + 14vN −1 ) = vN −2 − 2vN −1 + vN , 12 all formulae are O(h6 ), or written in matrix form: M V 00 = AV + H, where
A=
−2 1 0 .. . 0
M=
(2.5)
(2.6)
(2.7)
00 0 v1 .. v200 −2 1 . .. 00 .. .. .. and V = , . . . 0 . v 00 N −2 1 −2 1 00 vN −1 ··· 0 1 −2 (N −1)×(N −1) 14 −5 4 −1 0 · · · 0 .. 1 10 1 0 0 . 0 1 10 1 0 ··· 0 2 h .. .. .. , . . . 12 . .. 0 1 10 1 0 0 0 1 10 1 0 · · · 0 −1 4 −5 14 (N −1)×(N −1) v1 v0 v2 0 .. H = . , and V = ... . vN −2 0 vN −1 vN 1
0
···
Remark 2.1. Equation (2.5) and (2.6) are not compact finite difference formula, since at the left side they use four points. And the coefficients of them can be decided in a similar way we did for standard compact finite difference formula. Since the boundary values are known, we choose v0 = va and vN = vb in vector H. Write equation (2.4) into matrix form: −V 00 = F,
(2.8)
Compact Finite Difference Method for IDE
5
where F is a vector, its length is N − 1 and entries are f (xi ), i = 1, · · · , N − 1. By equation (2.7) and (2.8), we get −M V 00 = M F = −(AV + H), or AV = −(M F + H),
(2.9)
Let b = −(M F + H) which is known, so we obtain AV = b, Note the band width of matrix A is into A = LU , where 1 0 0 ··· −1 1 0 2 . . .. .. .. L= . 0 . .. −(N −3) 1 N −2 −(N −2) 0 ··· 0 N −1
(2.10)
three, and matrix A can be factored exactly 0 −2 1 .. −3 0 . 2 .. , U = 0 . 0 . . 0 . 1 0 ···
0
···
1 .. .
..
0
−(N −1) N −2
.
0
0 .. . 0 1 −(N ) N −1
We can use the decomposition to solve the matrix equations AV = (LU )V = L(U V ) = b,
(2.11)
by first solving for the vector Y such that LY = b,
(2.12)
i.e. y 1 = b1 , y i = b i −
−(i − 1) yi−1 , i = 2, · · · , N − 1, i
(2.13)
and then solving UV = Y
(2.14)
i.e. vN −1 = where
−(N − 1) i y N , vi = (vi+1 − yi ), i = (N − 2), · · · , 1, N i+1
(2.15)
y1 Y = ... , yN −1
So equation (2.11) can be solved for V very efficiently by LU decomposition. For the implementation issue, we can pre-factor the compact finite difference matrix so that we can save time on decomposing matrix A into L and U . We see that errors decrease in the above recurrences by the fact that |Li,i−1 | < 1, i = 2, · · · , N − 1, and |Ui,i | > 1, i = 1, · · · , N − 1. From the above, we can see that compact finite difference method for differential equation (2.3) can obtain high order of accuracy, while it only uses the time complexity of N to solve the matrix equations.
6
Jichao Zhao and Robert M. Corless
3. Compact Finite Difference Method for Second Order Integro-Differential Equations In this section, we will show how to use high order compact finite difference method for second order integro-differential equations. 3.1. Integro-Differential Equations with Typical Boundary Conditions For the integro-differential equation (1.2) with the typical boundary conditions (1.3), we need numerical integration formula to approximate the integral part of integro-differential equations. Gauss rules of numerical integration involve choosing the location of the nodes as well as the weights in an effort to reduce the number of function evaluations required to achieve a particular precision for the result. Generally speaking, the Gauss n point rule is exact on polynomials of degree 2n − 1, and requires n function evaluations. We choose Gauss 3 point rule in this section, At the point xi we apply Gauss 3 point rule of numerical integration and the three knots are r r h 3 h 3 h −1 0 1 ) , ωi = xi + , ωi = xi + (1 + ) , (3.1) ωi = xi + (1 − 5 2 2 5 2 ωi−1 , ωi0 and ωi1 ∈ [xi , xi+1 ], i = 0, · · · , N − 1, then we can have Z b K(xi , s)v(s)ds ≈ Qhi (Kv) =
N −1 h X 5 ( K(xi , ωi−1 )Pi3 (v, ωi−1 ) 2 i=0 9
a
(3.2)
8 5 K(xi , ωi0 )Pi3 (v, ωi0 ) + K(xi , ωi1 )Pi3 (v, ωi1 ), 9 9 where we replace v(x) by its interpolating cubic polynomial Pi3 (v, x) constructed using the four point values (xi−1 , vi−1 ), (xi , vi ), (xi+1 , vi+1 ), (xi+2 , vi+2 ) on subinterval [xi , xi+1 ], i = 1, · · · , (N − 2) [6]. For i = 0 and i = N − 1, we can choose P03 (v, x) = P13 (v, x) and PN3 −1 (v, x) = PN3 −2 (v, x) for subinterval [x0 , x1 ] and [xN −1 , xN ] respectively. And the global truncation error ei ≈ O(h4 ). Combining (1.2), (2.7) and (3.2), we get the following matrix equation +
M V 00 = AV + H = M (Qh (KV ) − F (V )), where
(3.3)
Qh1 f (x1 , v1 ) .. Qh = ... and F (V ) = . . h f (xN −1 , vN −1 ) QN −1
Let G(V ) = M (Qh (KV ) − F (V )) − H, so equation (3.3) can be written as AV = G(V ).
(3.4)
Compact Finite Difference Method for IDE
7
Remark 3.1. The global truncation error for the equation (3.4) is O(h6 ) since 2 matrix M itself has the common coefficient h12 . 3.2. Algorithm We solve the equation (3.4) by the following iterations [10] AV m+1 = G(V m ), m = 0, 1, · · · ,
(3.5)
where G(V m ) = M (Qh (KV m ) − F (V m )) − H. (3.6) m m+1 Suppose that we know the value V and want to compute V . Before we start to compute, we need to give initial values for V , say V 0 . We use the formula (3.2) and (3.3) to get a vector Qh (KV m ) − F (V m ), use the notation QF m to stand for it. Multiply QF m by the matrix M and subtract vector H, then we can obtain G(V m ) and the cost is 3(N − 1) + 2 of multiplications since matrix M is a sparse matrix. Finally, we solve equation (3.5) for V m+1 by (2.11∼2.15) in section 2 and the cost is only 2(N − 1). So the total cost is 5N − 3 of multiplications. In [6], the authors use finite difference scheme to solve integro-differential equation (1.2) with fourth order for interior points, but they need to factorize a pentagonal matrix into two tridiagonal matrices equations after they discretize the integrodifferential equation, then they transform them into upper diagonal equations by Gauss elimination. Their method uses more operations and more complicated than ours. Moreover, our compact finite difference method can obtain error O(h5 ) for the whole interval [a, b], while in [6] they only obtain the accuracy of O(h2 ) at the points x1 and xN −1 and O(h4 ) for interior mesh points. 3.3. Error Estimates To simplify, we use R = [a, b]. We define the discrete norms: k W k2h,0 =
N −1 X
hwi2 , | W |2h,1 =
i=1
N X (wi − wi−1 )2 i=1
h
,
for ∀ W = (w0 , w1 , · · · , wN )T ∈ RN +1 . Let R0N +1 = {∀ W ∈ RN +1 : w0 = wN = 0}, which is a subspace of RN +1 . It is not hard to know that ˜ ,W ˜ ), ∀ W ∈ RN +1 , k W k2h,0 = (hW (3.7) −1 ˜ ˜ | W |2h,1 = ( AW , W ), ∀ W ∈ R0N +1 , (3.8) h ˜ = (w1 , · · · , wN −1 )T . Define the continhere (·, ·) stands for inner product, and W 2 1 uous L and H norms as usual: Z b Z b k w k2 = w2 (x)dx, k w k21 = (wx2 (x) + w2 (x))dx, a
a 1
We define an interpolant operator Πh : H (R) → Rh , and Rh is a polynomial space. It is obvious that w(xi ) = Πh w(xi ), i = 0, 1, · · · , N , for ∀ w(x) ∈ H 1 (R).
8
Jichao Zhao and Robert M. Corless
By the approximation theory, in general we can obtain the following well known estimate for the interpolant: k w − Πh w k + h k w − Πh w k1 ≤ Chr k w kr , r > 1,
(3.9)
here the r is the degree of the approximate polynomials. By compact finite difference method, we obtain a series of discrete points vi to approximate the true solution ui at those partition points xi , i = 1, · · · , N −1, We can construct the interpolating polynomial Qh (x) by these discrete points (xi , vi ). Let eh = Πh u(x)−Qh (x), so ei = eh (xi ) = Πh u(xi ) − Qh (xi ) = ui − vi , i = 1, · · · , N − 1. Write them in the matrix form: x0 ua ua 0 x1 u1 v1 e1 .. .. .. X = . , U = . , V = . , and E = ... . xN −1 uN −1 vN −1 eN −1 xN ub ub , 0 For the discrete norms and continuous norms, we can easily obtain the following theorem (refer to [5] and [9]): Theorem 3.2. The norms k · k2h,0 and | · |2h,1 are equivalent to k · k2 and k · k21 respectively, namely, there exists positive constants C1 , C2 , C3 , and C4 independent of h such that C1 k W k2h,0 ≤ k Πh w k2 ≤ C1 k W k2h,0 , C3 | W
|2h,1
≤ k Πh w
k21
≤ C4 | W
|2h,1 ,
(3.10) (3.11)
for W = (w0 , w1 , · · · , wN )T ∈ RN +1 and w(x) ∈ H 1 (R). By the equivalent relations (3.10), (3.11) and inequality k w k2 ≤ k w k21 in Sobolev spaces, we can obtain: C k W k2h,0 ≤ | W |2h,1 ,
(3.12)
here C is a positive constant. We use a similar way with [6] to derive the error estimates for compact finite difference method of integro-differential equations, to this end, we first give another two assumptions which are similar to assumptions 1 and 2 in section 1: (V ) Assumption 3: The forcing function F (V ) and its partial derivatives ∂F∂u ˜ ,W ˜ ) is nonpositive for ∀ W ˜. are continuous and inner product ( h1 M ∂F∂u(η) W
Assumption 4: The kernel function K(X, S) satisfies the following positive definite property: Z b 1 ˜ ) > 0. ( M K(X, S)w(S)dS, W (3.13) h a
Compact Finite Difference Method for IDE
9
The assumptions 3 and 4 are inner product forms of assumptions 1 and 2 at partition points if the matrix M is an identity matrix. Theorem 3.3 (Error estimates). If the assumptions 3 and 4 are satisfied and the kernel functions K(x, s) is smooth enough, then the approximate solution Qh (x) converges to the solution u(x), and holds the inequality k u − Qh k ≤ Ch5 ,
(3.14)
where the C is any postive constant. Proof. Write u − Qh = u − Πh u + Πh u − Qh . By (3.9), we have k u − Πh u k ≤ Ch5 k u k5 ,
(3.15)
then we only need to estimate for eh = Πh u − Qh . Note that (Πh u − Qh )|x=xi = ui − vi = ei at partition point xi , i = 1, · · · , N − 1. Use ei + vi to replace u(x) at point xi in equation (1.2), and rewrite it into: Z b d2 ei − 2 + K(xi , s)eh (s)ds = f (xi , ei + vi ) − f (xi , vi ) + Φ(xi , vi ), (3.16) dx a where Z b d2 vi − Φ(xi , vi ) = f (xi , vi ) + K(xi , s)Qh (s)ds, (3.17) dx2 a and matrix form of equation (3.16) is: Z b d2 E K(X, S)eh (S)dS = F (E + V ) − F (V ) + Φ(V ). (3.18) − 2 + dx a By mean-value theorem, we have ∂F (η) E, (3.19) ∂u where η = (η1 , · · · , ηN −1 )T . Multiply both sides of equation (3.18) by matrix h1 M , and use the equation (2.7) and (3.8), then we obtain Z b 1 1 ∂F (ξ) 1 1 K(X, S)eh (S)dS, E) = ( M E, E)+( M Φ(V ), E), −( AE, E)+ ( M h h h ∂u h a (3.20) Z b 1 1 1 ∂F (ξ) | E |2h,1 + ( M E, E) + ( M Φ(V ), E), K(X, S)eh (S)dS, E) = ( M h h ∂u h a (3.21) By assumption 3, assumption 4, inequality (3.12) and Cauchy’s theorem, we have F (E + V ) − F (V ) =
1 1 C 1 k M Φ k2 + k E k2h,0 . C k E k2h,0 ≤ ( M Φ(V ), E) ≤ h 2C h 2 Using theorem 3.2, then we obtain 1 1 k eh k2 ≤ 2 k M Φ k2 , C h
(3.22)
(3.23)
10
Jichao Zhao and Robert M. Corless
Next, we will estimate for the k h1 M Φ k. By equations (3.2), (3.3) and (3.17), we obtain Z b 00 | M Φ(V ) | ≤ | M F (V ) + M V − M K(X, S)P3 (V, S)dS | Z
a b
+ |
M K(X, S)(P3 (V, S) − Qh (S))dS | a 6
≤ Ch
(3.24)
Combining inequality (3.23), (3.24) and (3.15), we can get inequality (3.14). This completes the proof. ¤ 3.4. Experiments We use the powerful software Maple to implement our algorithm. To measure the accuracy of our compact finite difference method, we use the residual error and exact error. The residual error R(xi ) is defined by Z b d2 vi R(xi ) = − 2 + K(xi , s)v ds − f (xi , vi ), 1 ≤ i ≤ N − 1, dx a and we use the higher order Taylor expansion (we use seventh order) to approxiRb 2 mate ddxv2i and Gauss 4 point rule of numerical integration formula for a K(xi , s)v ds in the above definition. If we know the exact solution for it, we can define the exact error as the absolute difference between the true value and the approximation value at every partition point, i.e. E(xi ) = | u(xi ) − v(xi ) |, 1 ≤ i ≤ N − 1. Example 1: d2 u(x) − + dx2
Z
1
se(−s) cos (x)u(s) ds = −2u(x) − x, −1 ≤ x ≤ 1,
(3.25)
−1
u(−1) = −1, u(1) = 1, the kernel function K(x, s) = cos(x)se(−s) , f (x, u) = −2u(x) − x, and the exact solution is: √
√
u(x) = C1 e 2x + C2 e− 2x + C3 cos (x) − .5x, C1 = .4206 · · · , C2 = −.3545 · · · , C3 = −.2662 · · · . We obtain the Fig. 1 for the numerical solutions of our high accurate compact finite difference method and this figure is for the fixed spacial step h. From the Fig. 1, we can see the exact error is between O(h7 ) and O(h9 ), while residual error is between O(h5 ) and O(h8 ). Next, we will show how the exact and residual error change with varying h by giving the Fig. 2. From the Fig. 2 and following numerical experiments, it seems that residual error overestimates the error and the numerical solutions by compact finite difference method is O(h5 ).
Compact Finite Difference Method for IDE
11
1e-06
Error
1e-07
1e-08
-0.8
-0.4
0
0.4
0.8
x Exact error Residual Error
Figure 1. Exact and residual err. for example 1, h = 0.1
To compare with the finite difference method used in [6], we use compact finite difference method for the same nonlinear integro-differential equation. Example 2: −
d2 u(x) + dx2
Z
1
K(x, s)u(s)ds = e−u(x) , 0 ≤ x ≤ 1
(3.26)
0
where the kernel K(x, s) is given by ½ x(1 − s), when 0 ≤ x ≤ s, K(x, s) = s(1 − x), when s ≤ x ≤ 1, and the boundary conditions of equation (3.26) are as follows: u(0) = 0, u(1) = 1.
(3.27)
12
Jichao Zhao and Robert M. Corless
.1e-3
1e-05 1e-06
1e-07 Error 1e-08
1e-09
1e-10 1e-11 .5e-2
.1e-1
.5e-1
.1
h Exact error Residual error
Figure 2. Exact and residual error with varying h for example 1
1 Let h = 10 , set initial value vi0 = 0, i = 1, · · · , 9. After eight iterations, we get the following numerical solutions :
xi vi Ri (10−6 ) xi vi Ri (10−6 )
0.1 0.128299 9.81633 0.5 0.566007 1.71506
0.2 0.247973 6.14125 0.6 0.66109 1.26291
0.3 0.360189 3.22083 0.7 0.751704 0.75500
0.4 0.465928 2.45317 0.8 0.838242 0.19599
0.5 0.566007 1.71506 0.9 0.920964 8.78495
Table 3.1 The above table also compact finite difference method can even obtain O(h6 ) on all the partition points, while finite difference method used in [6] only obtain O(h2 ) for the two points x1 and x9 . Still set vi0 = 0, i = 1, · · · , N − 1 and let spacial step h varies, we get the following table :
Compact Finite Difference Method for IDE Step h 0.1 0.05 0.025 0.0125 0.01 0.005 0.0025 0.001
Iteration numbers m 8 10 12 14 15 17 19 22
13
Maximum absolute residual errors 9.81633e − 06 4.86653e − 07 3.55414e − 08 2.50501e − 09 1.04889e − 09 6.77362e − 11 4.28567e − 12 1.10441e − 13
Table 3.2 From the above table, it shows that compact finite difference method can obtain error O(h5 ) at all the mesh points in the interval [a, b], but when h is too small, rounding errors take over. Last we will try different initial values vi0 , i = 1, · · · , N − 1, and see whether 1 , and we get the following table: our method depends on them or not. Let h = 50 tests test 1 test 2 test 3 test 4 test 5 test 6
Initial values vi0 vi0 = i/50, i = 1, · · · , 49 vi0 = 1 − i/50, i = 1, · · · , 49 vi0 = 0, i = 1, · · · , 49 vi0 = 100, i = 1, · · · , 49 vi0 = −8, i = 1, · · · , 49 vi0 = −10, i = 1, · · · , 49
Iteration numbers m 12 12 13 14 16 >100
Max absolute residual err. 1.50966e − 08 1.50966e − 08 1.50966e − 08 1.50966e − 08 1.50966e − 08 NaN
Table 3.3 It seems that the compact finite difference method converges very fast for test 1 ∼ test 5 and doesn’t depend too much on initial values. But for test 6, it doesn’t converge, the reason is that the term e−v(x) at the right side of the equation (1.2) is a large number when V 0 = −10. Suitably choosing V 0 can make us get desirable solutions V in a shorter time, i.e., smaller iterations m and guarantee this method converges. From the above numerical experiments, we can easily see that the compact finite difference method indeed obtain O(h5 ) and this method also converges fast. 3.5. Dealing with Different Boundary Conditions In this part, suppose we still have integro-differential equation (1.2), but with different boundary conditions instead of typical boundary conditions (1.3). When we have the following boundary conditions: u(a) = ua , u0 (b) = udb ,
(3.28)
u0 (a) = uda , u(b) = ub ,
(3.29)
or
we can deal with them almost in the same way as we did in section 2.2. The only change we need to make is that to use the following sixth order equation h2 127 00 86 00 257 00 461 00 0 (− vN −4 + vN vN −2 + v ) = vN −2 − vN −1 + hvN , (3.30) −3 − 12 30 5 10 15 N −1
14
Jichao Zhao and Robert M. Corless
for the differential boundary condition of (3.28) or h2 127 00 86 00 257 00 461 00 (− v + v3 − v + v ) = −v1 + v2 − hv00 , (3.31) 12 30 4 5 10 2 15 1 for the differential boundary condition of (3.29), instead of equation (2.6) and equation (2.5) respectively. Matrix M , A and vector H need to be modified accordingly and the algorithm is same with typical boundary conditions. In some conditions, we may meet with the following boundary conditions u(a) = ua , u0 (a) = uda ,
(3.32)
For this case, the boundary conditions are quite different with the boundary conditions we have mentioned before, since we only know the informations at one of two end points. The equation (2.6) can not be used since we don’t know the value of vN . We need one more equation except equation (2.5) and equations (2.1), so we try to use the known values, u0 (a) and u(a) together, to get another equation: h2 97 00 132 00 207 00 502 00 ( v − v + v − v ) = −2v1 + 2v0 + 2hv00 , (3.33) 12 15 4 5 3 5 2 15 1 also the above equation can obtain O(h6 ) order. Same thing for the following boundary conditions: u(b) = ub , u0 (b) = udb ,
(3.34)
We use the following equation except equation (2.1) and equation (2.6): h2 97 00 132 00 207 00 502 00 0 ( v − v + v − v ) = −2vN −1 + 2vN − 2hvN . (3.35) 12 15 N −4 5 N −3 5 N −2 15 N −1 For the following boundary conditions: u0 (a) = uda , u0 (b) = udb ,
(3.36)
we meet some problems. We may think of using equations similar to (3.31) and (3.30) to deal with boundary conditions (3.36), but no matter what accuracy of the formula will be, when we write them with equations (2.1), i = 1, · · · , N − 1 together into matrix form, compact finite difference matrix A always be: −1 1 0 ··· 0 0 .. .. 1 −2 1 . . . . . .. .. .. 0 0 0 A= , . .. .. .. .. . . . 0 0 1 −2 1 0 0 ··· 0 1 −1 but note det(A) = 0. So it means for boundary conditions (3.36), we need more conditions to solve integro-differential equation (1.2). In [1], it mentions the following general boundary conditions: a1 u(a) + b1 u(b) = λ1 , a2 u0 (a) + b2 u0 (b) = λ2 ,
(3.37)
Compact Finite Difference Method for IDE
15
where a1 , a2 , b1 , b2 , λ1 , and λ2 are known constants. Using Taylor expansion, we can use a linear combination of u(a), u1 , u2 , u3 , u4 and u5 to replace u0 (a) in (3.37) with sixth order of accuracy. Same way to replace u0 (b) by a linear combination of u(b), uN −1 , uN −2 , uN −3 , uN −4 and uN −5 . Solving for u(a) and u(b), we can get: u(a) = φ1 (u1 , u2 , u3 , u4 , u5 , uN −5 , uN −4 , uN −3 , uN −2 , uN −1 ), u(b) = φ2 (u1 , u2 , u3 , u4 , u4 , uN −5 , uN −4 , uN −3 , uN −2 , uN −1 ). Since we have the above relations, we can just deal with them in the same way as the typical boundary conditions by compact finite difference method. At the end of this part, we will give a procedure to show how to handle the case when we have the conditions are not limited on boundary points. Take a simple example to explain, suppose we have the following conditions: u(a) = ua , u(c) = uc ,
(3.38)
where ua and uc are two already known constants, and c ∈ (a, b). The procedure to handle them are as follows: 1. Solve the integro-differential equation on the interval [a, c] by compact finite difference method since we know the boundary conditions (3.38). So we can obtain all the numerical values for the mesh points in the subinterval [a, c]. 2. Construct another boundary condition u0 (c) by the linear combination of five mesh points in interval [a, c] closer to point c with the accuracy of sixth order by Taylor expansion. 3. Use the values of u0 (c) and u(c) as the boundary conditions to solve the integro-differential equation on the interval [c, b] by compact finite difference method. Then we obtain all the numerical values for the whole interval [a, b] with the condition (3.38). We can solve for other types of mixed boundary and interior points in the similar way.
4. Solving First Order Integro-Differential Equations We consider how to use compact finite difference method to solve the following first order integro−differential equation: Z b du(x) + K(x, s)u(s)ds = f (x, u), a ≤ x ≤ b, (4.1) − dx a with boundary condition u(a) = u0 .
(4.2)
The standard compact finite difference formula for first derivatives is: 0 0 h(vi−1 + 4vi0 + vi+1 ) = 3(−vi−1 + vi+1 ),
(4.3)
16
Jichao Zhao and Robert M. Corless
and the error in this formula is O(h5 ). To consistent with the truncation error of O(h6 ) for second derivatives in section 3 and get a better property of matrix A, we use the following formula for first derivatives: h (−23v10 + 39v20 − 21v30 + 5v40 ) = v0 − 2v1 + v2 , (4.4) 12 9 1 h (− v10 − 5v20 + 12v30 − 3v40 + v50 ) = v1 − 2v2 + v2 , (4.5) 12 2 2 h 1 0 1 0 0 0 ( v − 7vi−1 + 7vi+1 − vi+2 ) = vi−1 − 2vi + vi+1 , i = 3, · · · , N − 2, (4.6) 12 2 i−2 2 h 1 0 9 0 0 0 0 (− vN −4 + 3vN −3 − 12vN −2 + 5vN −1 + vN ) = vN −2 − 2vN −1 + vN , (4.7) 12 2 2 h 1 0 5 0 19 0 9 0 (− v + v − vN −1 − vN ) = vN −1 − vN , (4.8) 12 2 N −3 2 N −2 2 2 6 all formulae are O(h ), and can be obtained by Taylor expansion. Or rewrite them in matrix form: M V 0 = AV + H, (4.9) where A=
−2
1
1
−2 .. .
0 .. . 0
···
−23 9 − 2 0 . h .. M= 12 0 . .. ··· 0
H
0 0 v1 .. v20 1 . .. .. and V 0 = ... , . . 0 v 0 N −1 1 −2 1 0 vN 0 1 −1 N ×N 39 −21 5 0 ··· 0 .. −5 12 −3 − 21 0 . 1 1 −7 0 7 −2 0 2 .. .. .. .. .. .. . . . . . . , 1 1 −7 0 7 −2 0 2 1 0 −7 0 7 − 21 2 9 0 − 12 3 −12 5 2 ··· 0 −5 21 −39 23 N ×N v1 v0 v2 0 .. = . and V = ... . vN −1 0 vN 0 0
···
Compact Finite Difference Method for IDE Matrix A can 1 −1 2 L= 0 .. . 0
be factored exactly into LU : 0 0 ··· 0 −2 1 .. −3 0 1 0 . 2 .. .. .. . .. . . . 0 , U = 0 . −(N −2) .. 1 0 N −1 −(N −1) 0 ··· ··· 0 1 N
17
0
···
1 .. .
..
0
−N N −1
.
0
0 .. . . 0 1 −1 N
For the integral part of (4.1), we use the same way as we did for second order integro-differential equation. Then we can solve the first order integro−differential problem (4.1) and (4.2) with high order of accuracy by compact finite difference method. Finally we give an example to show our compact finite difference method works for first order integro−differential equation. Example 3: Z 1 du(x) − + e(−s) sin xu(s)ds = u(x) − 2x, 0 ≤ x ≤ 1, (4.10) dx 0 u(0) = 0, The kernel function K(x, s) = sin(x)e(−s) , f (x, u) = u(x) − 2x, and the exact solution is: u(x) = C1 sin(x) + 2x − 2 + C2 cos(x) + C3 e(−x) , C1 = .0686 · · · , C2 = −.0686 · · · , C3 = 2.0686 · · · , By compact finite difference method, we get the Fig. 3 and 4 for example 3.
5. Solving the System of Integro-Differential Equations When we have a system of integro-differential equations, we can also use compact finite differene method to solve them. To make life easier, we take J, the number of integro-differential equations, as the two and both of them are second order integro-differential equations, to explain how to solve the system. In the similar way, we can solve for the case J > 2 or the system has first order integro-differential equations. In this section, we will focus on the following system of two second order integro-differential equations: Z b Z b d2 u1 (x) − + K11 (x, s)u1 (s)ds + K12 (x, s)u2 (s)ds = f1 (x, u1 , u2 ), (5.1) dx a a Z Z b b d2 u2 (x) + K21 (x, s)u1 (s)ds + K22 (x, s)u2 (s)ds = f2 (x, u1 , u2 ), (5.2) − dx a a a ≤ x ≤ b,
18
Jichao Zhao and Robert M. Corless
7e-07
4e-07
Error 2e-07
1e-07 0.2
0.4
0.6
0.8
x Exact error Residual error
Figure 3. Exact and residual err. for example 3, h = .1 with boundary conditions u1 (a) = u1a , u1 (b) = u1b , u2 (a) = u2a , u2 (a) = u2b ,
(5.3) (5.4)
And the algorithm to solve the above (5.1∼5.4) is given below: 1. By compact finite difference method, we can write the differential parts of equations (5.1) and (5.2) into matrix form: M V100 = AV1 + H1 , M V200 = AV2 + H2 , where the elements of vector V1 and V2 correspond to partition points of functions u1 (x) and u2 (x), please refer to equation (2.7) of section 2 for details. 2. We can obtain the two following matrix equations after we combine the above two matrix equations and Gauss 3 point rule for the integral parts of equations (5.1) and (5.2): AV1 = G(V1 , V2 ), AV1 = G(V1 , V2 ).
(5.5)
Compact Finite Difference Method for IDE
19
1e-06
1e-07
1e-08
Error
1e-09
1e-10
1e-11
.5e-2
.1e-1
.2e-1
.5e-1
.1
h Exact error Residual error
Figure 4. Exact and residual err. with varying h for example 3
3. Set initial values for vector V10 and V20 , then compute new V1m+1 and V2m+1 using already known values of V1m and V2m by equations (5.5), m = 0, 1, · · · , until the difference, both V1m+1 − V1m and V2m+1 − V2m , is smaller than a given error. The following example, a system of integro-differential equations, is solved by compact finite difference method. Example 4: d2 u1 (x) − + dx2 −
Z
Z
1
sin (x)su1 (s)ds + 2 0
d2 u2 (x) + dx2
Z
Z
1
su1 (s)ds + 0
1
e−s u2 (s)ds = −u1 (x) − e−x ,
0 1
cos (x)e−s u2 (s)ds = −u2 (x) − 4,
0
0 ≤ x ≤ 1,
20
Jichao Zhao and Robert M. Corless
with boundary conditions u1 (0) = 0, u1 (1) = 1, u2 (0) = 0, u2 (1) = 1,
(5.6) (5.7)
and the analytic solution is: 1 u1 (x) = C10 + C11 sin(x) + C12 e(x) + C13 e(−x) − xe(−x) , 2 u2 (x) = C20 + C21 cos(x) + C22 e(x) + C23 e(−x) , C10 = −.0763 · · · , C11 = −.1410 · · · , C12 = .5747 · · · , C13 = −.4983 · · · , C20 = −4.2820 · · · , C21 = −.0190 · · · , C22 = 1.5784 · · · , C23 = 2.7226 · · · . By compact finite difference method, we get the numerical solutions for u1 and u2 , please see Fig. 5 ∼ 8 for them.
6. Conclusions Compact finite difference method can solve very general integro-differential equations with high accuracy, while only use time complexity O(N ) for solving the matrix equations. It is easy to extend our work for Volterra integro-differential equation.
References [1] A. Karameta, and M. Sezer, A Taylor collocation method for the solution of linear integro-differential equations, Intern. J. Computer Math., 2002, Vol. 79(9), pp. 9871000. [2] E. Liz, and J.J. Nieto, Boundary value problems for second order integro−differential equations of Fredholm type, J. Comput. Appl. Math. 72 (1996) 215−225. [3] J. Zhao, and T. Zhang, A highly accurate derivative recovery formula to integrodifferential equations, Numerical Mathematics a Journal of Chinese Universities, 2004 Vol.26(1), pp.81-90. [4] H. Sun, and J. Zhang, A high order compact boundary value method for solving one dimensional heat equations, Technical Report No. 333-02, University of Kentucky, 2002. [5] H.Yang, Numerical analysis of American options, Ph.D. Thesis, University of Alberta, 2001. [6] K. Sty´s, and T. Sty´s, A higher-order finite difference method for solving a system of integro-differential equations, Journal of Computational and Applied Mathematics, 126(2000) 33-46. [7] M.O. Ahmed, An exploration of compact finite difference methods for the numerical solution of PDE, Ph.D. thesis, the University of Western Ontario, 1997.
Compact Finite Difference Method for IDE
21
1e-06
8e-07
6e-07 Error
4e-07
3e-07 0.2
0.4
0.6
0.8
x Exact error Residual error
Figure 5. Exact and residual err. for u1 (x), h = 0.1
[8] R. Cont, P. Tankov and E. Voltchkova, Option pricing models with jumps: integrodifferential equations and inverse problems. European Congress on Computational Methods in Applied Sciences and Engineering, 24-28 July (2004) 1-20. [9] R.H. Li, Z.Y. Chen and W. Wu, Generalized difference methods for differential equations (Numerical analysis of finite volume methods), Marcel Dekker, Inc., New York, 2000. [10] R.S. Varga, Matrix iterative analysis, Berlin; New York: Springer, c2000. [11] R.M. Corless, and J. Rokicki, The symbolic generation of finite-difference formulas. Zeitschrift fur Angewandte Mathematik und Mechanik (ZAMM) Proceedings. ICIAM/GAMM 95, Numerical Analysis, Scientific Computing, Computer Science. Hamburg, Germany: 3 July 1995-7 July 1995. G. Alefeld, O. Mahrenholtz, and R. Mennicken. ZAMM, Vol. 76, Supplement. Germany: Akademie Verlag. 381-382.
22
Jichao Zhao and Robert M. Corless
2e-06
1e-06
Error 7e-07
4e-07
0.2
0.4
0.6
0.8
x Exact error Residual error
Figure 6. Exact and residual err. for u2 (x), h = 0.1
[12] R.M. Corless, J. Rokicki and J. Zhao, FINDIF, a routine for generation of finite difference formulae, Share Library package 1994, and upgraded to n dimensions for “iguana”. [13] R. Natalini, M. Briani and C.L. Chioma, Finite differences schemes for integrodifferential equations in financial mathematics, IAC-CNR & University, 2003. [14] T. Zhang, Finite element methods for evolution integro-differential equations, Northeastern University Press, Shenyang of China, 2001. [15] W. Liu, The exponential stabilization of the high dimensional linear system of thermoviscoelasticity. J. Math. Pures Appl., 77 (1998) 355-386.
Compact Finite Difference Method for IDE
23
1e-06
1e-07
1e-08 Error
1e-09
1e-10
1e-11 .5e-2
.1e-1
.2e-1 h
.5e-1 Exact error Residual error
Figure 7. Exact and residual err. for u1 (x) with varying h Acknowledgment Zhao is thankful to Prof. T. Zhang for remarks on a draft of this paper. Jichao Zhao Department of Applied Mathematics, University of Western Ontario, Middlesex College, London, Ontario, Canada N6A 5B7 e-mail: [email protected] Robert M. Corless Department of Applied Mathematics, University of Western Ontario, Middlesex College, London, Ontario, Canada N6A 5B7 e-mail: [email protected]
.1
24
Jichao Zhao and Robert M. Corless
1e-06
1e-07 1e-08 Error 1e-09
1e-10 1e-11
.5e-2
.1e-1
.2e-1 h
.5e-1 Exact error Residual error
Figure 8. Exact and residual err. for u2 (x) with varying h
.1
1. Introduction Integro-differential equation (IDE) is an equation that the unknown function appears under the sign of integration and it also contains the derivatives of the unknown function. It can be classified into Fredholm equations and Volterra equations. The upper bound of the region for integral part of Volterra type is variable, while it is a fixed number for that of Fredholm type. For example, the following integro-differential equation belongs to Volterra type: Z x d2 u(x) + K(x, s)u(s)ds = f (x, u), a ≤ x ≤ b (1.1) − dx2 a and equation (1.2) is Fredholm integro-differential equation. In this paper we focus on Fredholm integro-differential equations, but all algorithms in our work can be applied to Volterra integro-differential equation with a little modification. If the kernel function holds K(x, s, u(s)) = K(x, s)u(s), we call the kernel is “split” type,
2
Jichao Zhao and Robert M. Corless
otherwise it is “unsplit” type. Our compact finite difference method works even for unsplit kernel, to simplify, we write the kernel function separately. Integro-differential equations arise from many branches of science, for example in control theory and financial mathematics (refer to [15], [8], [13], and [14]). Integro-differential equations are important, but they are hard to solve even numerically, so the progress on how to solve them is slow. In this paper, we mainly deal with the following integro-differential equation [6]: Z b d2 u(x) − + K(x, s)u(s)ds = f (x, u), a ≤ x ≤ b, (1.2) dx2 a where function f (x, u) can be a nonlinear function of u(x), and the typical boundary conditions: u(a) = ua , u(b) = ub ,
(1.3)
here ua and ub are already known boundary values. We assume the above equation (1.2) under the following two assumptions [6]: 1. The forcing function f (x, u) and all its partial derivatives are continuous , (x,u) and ∂f ∂u is nonpositive for a ≤ x ≤ b. 2. The kernel function K(x, s) satisfies the following positive definite property: Z bZ b K(x, s)ξ(x)ξ(s)dsdx > 0, (1.4) a
a
where ξ(x) is any continuous nonzero function, and holds Z bZ b | K(x, s) |2 dsdx < +∞. a
(1.5)
a
We know from the references [2], [6], and [14], the problem (1.2∼1.3) has a unique solution u(x) under the above two assumptions. In section 2, we give a brief introduction to a high accurate compact finite difference formula for second derivatives. In section 3, we show how to use compact finite difference method to solve the second order integro-differential equations with varying boundary conditions, also give the algorithm and error estimates for them. Numerical experiments in the last part of this section show that compact finite difference method can obtain fifth order of accuracy. In section 4, we give a high accurate compact finite difference method for first derivatives and use it to solve first order integro-differential equations, then followed by numerical examples to show this method also obtains fifth order of accuracy. As an expansion, we describe how to solve a system of integro-differential equations and an illustration is given in section 5. Conclusions are drawn in section 6. Note we use v(x) and u(x) to denote the numerical solution and exact solution of integro-differential equations respectively in the following sections.
Compact Finite Difference Method for IDE
3
2. Sixth Order Compact Finite Difference Formula Compact finite difference method is a special finite difference method which uses the values of the function and its derivatives only at three consecutive points. First discretize the interval [a, b], let spatial step h = 4x = |b−a| N , xi = a + ih, i = 0, 1, · · · , N − 1, N . 2.1. Standard Compact Finite Difference Formula The standard sixth order compact finite difference formula for second derivatives is d2 vi d2 vi+1 h2 d2 vi−1 ( + 10 + ) = vi−1 − 2vi + vi+1 , (2.1) 12 dx2 dx2 dx2 where vi = v(xi ), and the coefficients can be determined in the following way: 1. Write the desired compact finite difference formula with unknown coefficients: d2 vi−1 d2 vi d2 vi+1 + m0 2 + m1 ) = a−1 vi−1 + a0 vi + a1 vi+1 , (2.2) 2 dx dx dx2 here m−1 , m0 , m1 , a−1 , a0 , and a1 are the parameters to be decided.
h2 (m−1
2. Expand both sides of the equation (2.2) by Taylor series at the point xi with respect to the discretization parameter h, then collect them by the order of h. 3. We can obtain six equations by setting the coefficients of hj , j = 0, 1, · · · , 4, 5, equal zero. Finally solve the six equations for the six unknown parameters. From the above algorithm, we can see easily the accuracy is O(h6 ) for formula (2.1). For more details on how to generate compact finite difference formula, please refer to [7], [11], and [12]. Next we will show how to adjust compact finite difference formula for differential equations. 2.2. Solving Differential Equations by Compact Finite Difference Method Suppose we want to solve the following second order differential equation numerically for u(x): d2 u(x) = f (x), a ≤ x ≤ b, with u(a) = ua , u(b) = ub , (2.3) dx2 where f (x) is a known function. Write equation (2.3) in a discrete point form: −
d2 u(xi ) = f (xi ), 1 ≤ i ≤ N − 1, with u0 = ua , uN = ub . (2.4) dx2 Now the problem is that we know the values of second derivatives u00 (xi ) by equation (2.4) and need to solve for u(xi ), 1 ≤ i ≤ N − 1. So if we can find some relation between u00 (xi ) and u(xi ), then we can solve the equation (2.3) numerically. The tradition way is to use the second accurate central difference formula, while compact finite difference formula provides a higher accurate method to solve equation (2.3). We use standard sixth order compact finite difference formula (2.1), −
4
Jichao Zhao and Robert M. Corless
i = 2, · · · , (N − 2), for interior points. Since we have known the boundary conditions, we need to adjust compact finite difference formula for the boundary points. When i = 1, we use h2 (14v100 − 5v200 + 4v300 − v400 ) = v0 − 2v1 + v2 , 12 and when i = N − 1, h2 00 00 00 00 (−vN −4 + 4vN −3 − 5vN −2 + 14vN −1 ) = vN −2 − 2vN −1 + vN , 12 all formulae are O(h6 ), or written in matrix form: M V 00 = AV + H, where
A=
−2 1 0 .. . 0
M=
(2.5)
(2.6)
(2.7)
00 0 v1 .. v200 −2 1 . .. 00 .. .. .. and V = , . . . 0 . v 00 N −2 1 −2 1 00 vN −1 ··· 0 1 −2 (N −1)×(N −1) 14 −5 4 −1 0 · · · 0 .. 1 10 1 0 0 . 0 1 10 1 0 ··· 0 2 h .. .. .. , . . . 12 . .. 0 1 10 1 0 0 0 1 10 1 0 · · · 0 −1 4 −5 14 (N −1)×(N −1) v1 v0 v2 0 .. H = . , and V = ... . vN −2 0 vN −1 vN 1
0
···
Remark 2.1. Equation (2.5) and (2.6) are not compact finite difference formula, since at the left side they use four points. And the coefficients of them can be decided in a similar way we did for standard compact finite difference formula. Since the boundary values are known, we choose v0 = va and vN = vb in vector H. Write equation (2.4) into matrix form: −V 00 = F,
(2.8)
Compact Finite Difference Method for IDE
5
where F is a vector, its length is N − 1 and entries are f (xi ), i = 1, · · · , N − 1. By equation (2.7) and (2.8), we get −M V 00 = M F = −(AV + H), or AV = −(M F + H),
(2.9)
Let b = −(M F + H) which is known, so we obtain AV = b, Note the band width of matrix A is into A = LU , where 1 0 0 ··· −1 1 0 2 . . .. .. .. L= . 0 . .. −(N −3) 1 N −2 −(N −2) 0 ··· 0 N −1
(2.10)
three, and matrix A can be factored exactly 0 −2 1 .. −3 0 . 2 .. , U = 0 . 0 . . 0 . 1 0 ···
0
···
1 .. .
..
0
−(N −1) N −2
.
0
0 .. . 0 1 −(N ) N −1
We can use the decomposition to solve the matrix equations AV = (LU )V = L(U V ) = b,
(2.11)
by first solving for the vector Y such that LY = b,
(2.12)
i.e. y 1 = b1 , y i = b i −
−(i − 1) yi−1 , i = 2, · · · , N − 1, i
(2.13)
and then solving UV = Y
(2.14)
i.e. vN −1 = where
−(N − 1) i y N , vi = (vi+1 − yi ), i = (N − 2), · · · , 1, N i+1
(2.15)
y1 Y = ... , yN −1
So equation (2.11) can be solved for V very efficiently by LU decomposition. For the implementation issue, we can pre-factor the compact finite difference matrix so that we can save time on decomposing matrix A into L and U . We see that errors decrease in the above recurrences by the fact that |Li,i−1 | < 1, i = 2, · · · , N − 1, and |Ui,i | > 1, i = 1, · · · , N − 1. From the above, we can see that compact finite difference method for differential equation (2.3) can obtain high order of accuracy, while it only uses the time complexity of N to solve the matrix equations.
6
Jichao Zhao and Robert M. Corless
3. Compact Finite Difference Method for Second Order Integro-Differential Equations In this section, we will show how to use high order compact finite difference method for second order integro-differential equations. 3.1. Integro-Differential Equations with Typical Boundary Conditions For the integro-differential equation (1.2) with the typical boundary conditions (1.3), we need numerical integration formula to approximate the integral part of integro-differential equations. Gauss rules of numerical integration involve choosing the location of the nodes as well as the weights in an effort to reduce the number of function evaluations required to achieve a particular precision for the result. Generally speaking, the Gauss n point rule is exact on polynomials of degree 2n − 1, and requires n function evaluations. We choose Gauss 3 point rule in this section, At the point xi we apply Gauss 3 point rule of numerical integration and the three knots are r r h 3 h 3 h −1 0 1 ) , ωi = xi + , ωi = xi + (1 + ) , (3.1) ωi = xi + (1 − 5 2 2 5 2 ωi−1 , ωi0 and ωi1 ∈ [xi , xi+1 ], i = 0, · · · , N − 1, then we can have Z b K(xi , s)v(s)ds ≈ Qhi (Kv) =
N −1 h X 5 ( K(xi , ωi−1 )Pi3 (v, ωi−1 ) 2 i=0 9
a
(3.2)
8 5 K(xi , ωi0 )Pi3 (v, ωi0 ) + K(xi , ωi1 )Pi3 (v, ωi1 ), 9 9 where we replace v(x) by its interpolating cubic polynomial Pi3 (v, x) constructed using the four point values (xi−1 , vi−1 ), (xi , vi ), (xi+1 , vi+1 ), (xi+2 , vi+2 ) on subinterval [xi , xi+1 ], i = 1, · · · , (N − 2) [6]. For i = 0 and i = N − 1, we can choose P03 (v, x) = P13 (v, x) and PN3 −1 (v, x) = PN3 −2 (v, x) for subinterval [x0 , x1 ] and [xN −1 , xN ] respectively. And the global truncation error ei ≈ O(h4 ). Combining (1.2), (2.7) and (3.2), we get the following matrix equation +
M V 00 = AV + H = M (Qh (KV ) − F (V )), where
(3.3)
Qh1 f (x1 , v1 ) .. Qh = ... and F (V ) = . . h f (xN −1 , vN −1 ) QN −1
Let G(V ) = M (Qh (KV ) − F (V )) − H, so equation (3.3) can be written as AV = G(V ).
(3.4)
Compact Finite Difference Method for IDE
7
Remark 3.1. The global truncation error for the equation (3.4) is O(h6 ) since 2 matrix M itself has the common coefficient h12 . 3.2. Algorithm We solve the equation (3.4) by the following iterations [10] AV m+1 = G(V m ), m = 0, 1, · · · ,
(3.5)
where G(V m ) = M (Qh (KV m ) − F (V m )) − H. (3.6) m m+1 Suppose that we know the value V and want to compute V . Before we start to compute, we need to give initial values for V , say V 0 . We use the formula (3.2) and (3.3) to get a vector Qh (KV m ) − F (V m ), use the notation QF m to stand for it. Multiply QF m by the matrix M and subtract vector H, then we can obtain G(V m ) and the cost is 3(N − 1) + 2 of multiplications since matrix M is a sparse matrix. Finally, we solve equation (3.5) for V m+1 by (2.11∼2.15) in section 2 and the cost is only 2(N − 1). So the total cost is 5N − 3 of multiplications. In [6], the authors use finite difference scheme to solve integro-differential equation (1.2) with fourth order for interior points, but they need to factorize a pentagonal matrix into two tridiagonal matrices equations after they discretize the integrodifferential equation, then they transform them into upper diagonal equations by Gauss elimination. Their method uses more operations and more complicated than ours. Moreover, our compact finite difference method can obtain error O(h5 ) for the whole interval [a, b], while in [6] they only obtain the accuracy of O(h2 ) at the points x1 and xN −1 and O(h4 ) for interior mesh points. 3.3. Error Estimates To simplify, we use R = [a, b]. We define the discrete norms: k W k2h,0 =
N −1 X
hwi2 , | W |2h,1 =
i=1
N X (wi − wi−1 )2 i=1
h
,
for ∀ W = (w0 , w1 , · · · , wN )T ∈ RN +1 . Let R0N +1 = {∀ W ∈ RN +1 : w0 = wN = 0}, which is a subspace of RN +1 . It is not hard to know that ˜ ,W ˜ ), ∀ W ∈ RN +1 , k W k2h,0 = (hW (3.7) −1 ˜ ˜ | W |2h,1 = ( AW , W ), ∀ W ∈ R0N +1 , (3.8) h ˜ = (w1 , · · · , wN −1 )T . Define the continhere (·, ·) stands for inner product, and W 2 1 uous L and H norms as usual: Z b Z b k w k2 = w2 (x)dx, k w k21 = (wx2 (x) + w2 (x))dx, a
a 1
We define an interpolant operator Πh : H (R) → Rh , and Rh is a polynomial space. It is obvious that w(xi ) = Πh w(xi ), i = 0, 1, · · · , N , for ∀ w(x) ∈ H 1 (R).
8
Jichao Zhao and Robert M. Corless
By the approximation theory, in general we can obtain the following well known estimate for the interpolant: k w − Πh w k + h k w − Πh w k1 ≤ Chr k w kr , r > 1,
(3.9)
here the r is the degree of the approximate polynomials. By compact finite difference method, we obtain a series of discrete points vi to approximate the true solution ui at those partition points xi , i = 1, · · · , N −1, We can construct the interpolating polynomial Qh (x) by these discrete points (xi , vi ). Let eh = Πh u(x)−Qh (x), so ei = eh (xi ) = Πh u(xi ) − Qh (xi ) = ui − vi , i = 1, · · · , N − 1. Write them in the matrix form: x0 ua ua 0 x1 u1 v1 e1 .. .. .. X = . , U = . , V = . , and E = ... . xN −1 uN −1 vN −1 eN −1 xN ub ub , 0 For the discrete norms and continuous norms, we can easily obtain the following theorem (refer to [5] and [9]): Theorem 3.2. The norms k · k2h,0 and | · |2h,1 are equivalent to k · k2 and k · k21 respectively, namely, there exists positive constants C1 , C2 , C3 , and C4 independent of h such that C1 k W k2h,0 ≤ k Πh w k2 ≤ C1 k W k2h,0 , C3 | W
|2h,1
≤ k Πh w
k21
≤ C4 | W
|2h,1 ,
(3.10) (3.11)
for W = (w0 , w1 , · · · , wN )T ∈ RN +1 and w(x) ∈ H 1 (R). By the equivalent relations (3.10), (3.11) and inequality k w k2 ≤ k w k21 in Sobolev spaces, we can obtain: C k W k2h,0 ≤ | W |2h,1 ,
(3.12)
here C is a positive constant. We use a similar way with [6] to derive the error estimates for compact finite difference method of integro-differential equations, to this end, we first give another two assumptions which are similar to assumptions 1 and 2 in section 1: (V ) Assumption 3: The forcing function F (V ) and its partial derivatives ∂F∂u ˜ ,W ˜ ) is nonpositive for ∀ W ˜. are continuous and inner product ( h1 M ∂F∂u(η) W
Assumption 4: The kernel function K(X, S) satisfies the following positive definite property: Z b 1 ˜ ) > 0. ( M K(X, S)w(S)dS, W (3.13) h a
Compact Finite Difference Method for IDE
9
The assumptions 3 and 4 are inner product forms of assumptions 1 and 2 at partition points if the matrix M is an identity matrix. Theorem 3.3 (Error estimates). If the assumptions 3 and 4 are satisfied and the kernel functions K(x, s) is smooth enough, then the approximate solution Qh (x) converges to the solution u(x), and holds the inequality k u − Qh k ≤ Ch5 ,
(3.14)
where the C is any postive constant. Proof. Write u − Qh = u − Πh u + Πh u − Qh . By (3.9), we have k u − Πh u k ≤ Ch5 k u k5 ,
(3.15)
then we only need to estimate for eh = Πh u − Qh . Note that (Πh u − Qh )|x=xi = ui − vi = ei at partition point xi , i = 1, · · · , N − 1. Use ei + vi to replace u(x) at point xi in equation (1.2), and rewrite it into: Z b d2 ei − 2 + K(xi , s)eh (s)ds = f (xi , ei + vi ) − f (xi , vi ) + Φ(xi , vi ), (3.16) dx a where Z b d2 vi − Φ(xi , vi ) = f (xi , vi ) + K(xi , s)Qh (s)ds, (3.17) dx2 a and matrix form of equation (3.16) is: Z b d2 E K(X, S)eh (S)dS = F (E + V ) − F (V ) + Φ(V ). (3.18) − 2 + dx a By mean-value theorem, we have ∂F (η) E, (3.19) ∂u where η = (η1 , · · · , ηN −1 )T . Multiply both sides of equation (3.18) by matrix h1 M , and use the equation (2.7) and (3.8), then we obtain Z b 1 1 ∂F (ξ) 1 1 K(X, S)eh (S)dS, E) = ( M E, E)+( M Φ(V ), E), −( AE, E)+ ( M h h h ∂u h a (3.20) Z b 1 1 1 ∂F (ξ) | E |2h,1 + ( M E, E) + ( M Φ(V ), E), K(X, S)eh (S)dS, E) = ( M h h ∂u h a (3.21) By assumption 3, assumption 4, inequality (3.12) and Cauchy’s theorem, we have F (E + V ) − F (V ) =
1 1 C 1 k M Φ k2 + k E k2h,0 . C k E k2h,0 ≤ ( M Φ(V ), E) ≤ h 2C h 2 Using theorem 3.2, then we obtain 1 1 k eh k2 ≤ 2 k M Φ k2 , C h
(3.22)
(3.23)
10
Jichao Zhao and Robert M. Corless
Next, we will estimate for the k h1 M Φ k. By equations (3.2), (3.3) and (3.17), we obtain Z b 00 | M Φ(V ) | ≤ | M F (V ) + M V − M K(X, S)P3 (V, S)dS | Z
a b
+ |
M K(X, S)(P3 (V, S) − Qh (S))dS | a 6
≤ Ch
(3.24)
Combining inequality (3.23), (3.24) and (3.15), we can get inequality (3.14). This completes the proof. ¤ 3.4. Experiments We use the powerful software Maple to implement our algorithm. To measure the accuracy of our compact finite difference method, we use the residual error and exact error. The residual error R(xi ) is defined by Z b d2 vi R(xi ) = − 2 + K(xi , s)v ds − f (xi , vi ), 1 ≤ i ≤ N − 1, dx a and we use the higher order Taylor expansion (we use seventh order) to approxiRb 2 mate ddxv2i and Gauss 4 point rule of numerical integration formula for a K(xi , s)v ds in the above definition. If we know the exact solution for it, we can define the exact error as the absolute difference between the true value and the approximation value at every partition point, i.e. E(xi ) = | u(xi ) − v(xi ) |, 1 ≤ i ≤ N − 1. Example 1: d2 u(x) − + dx2
Z
1
se(−s) cos (x)u(s) ds = −2u(x) − x, −1 ≤ x ≤ 1,
(3.25)
−1
u(−1) = −1, u(1) = 1, the kernel function K(x, s) = cos(x)se(−s) , f (x, u) = −2u(x) − x, and the exact solution is: √
√
u(x) = C1 e 2x + C2 e− 2x + C3 cos (x) − .5x, C1 = .4206 · · · , C2 = −.3545 · · · , C3 = −.2662 · · · . We obtain the Fig. 1 for the numerical solutions of our high accurate compact finite difference method and this figure is for the fixed spacial step h. From the Fig. 1, we can see the exact error is between O(h7 ) and O(h9 ), while residual error is between O(h5 ) and O(h8 ). Next, we will show how the exact and residual error change with varying h by giving the Fig. 2. From the Fig. 2 and following numerical experiments, it seems that residual error overestimates the error and the numerical solutions by compact finite difference method is O(h5 ).
Compact Finite Difference Method for IDE
11
1e-06
Error
1e-07
1e-08
-0.8
-0.4
0
0.4
0.8
x Exact error Residual Error
Figure 1. Exact and residual err. for example 1, h = 0.1
To compare with the finite difference method used in [6], we use compact finite difference method for the same nonlinear integro-differential equation. Example 2: −
d2 u(x) + dx2
Z
1
K(x, s)u(s)ds = e−u(x) , 0 ≤ x ≤ 1
(3.26)
0
where the kernel K(x, s) is given by ½ x(1 − s), when 0 ≤ x ≤ s, K(x, s) = s(1 − x), when s ≤ x ≤ 1, and the boundary conditions of equation (3.26) are as follows: u(0) = 0, u(1) = 1.
(3.27)
12
Jichao Zhao and Robert M. Corless
.1e-3
1e-05 1e-06
1e-07 Error 1e-08
1e-09
1e-10 1e-11 .5e-2
.1e-1
.5e-1
.1
h Exact error Residual error
Figure 2. Exact and residual error with varying h for example 1
1 Let h = 10 , set initial value vi0 = 0, i = 1, · · · , 9. After eight iterations, we get the following numerical solutions :
xi vi Ri (10−6 ) xi vi Ri (10−6 )
0.1 0.128299 9.81633 0.5 0.566007 1.71506
0.2 0.247973 6.14125 0.6 0.66109 1.26291
0.3 0.360189 3.22083 0.7 0.751704 0.75500
0.4 0.465928 2.45317 0.8 0.838242 0.19599
0.5 0.566007 1.71506 0.9 0.920964 8.78495
Table 3.1 The above table also compact finite difference method can even obtain O(h6 ) on all the partition points, while finite difference method used in [6] only obtain O(h2 ) for the two points x1 and x9 . Still set vi0 = 0, i = 1, · · · , N − 1 and let spacial step h varies, we get the following table :
Compact Finite Difference Method for IDE Step h 0.1 0.05 0.025 0.0125 0.01 0.005 0.0025 0.001
Iteration numbers m 8 10 12 14 15 17 19 22
13
Maximum absolute residual errors 9.81633e − 06 4.86653e − 07 3.55414e − 08 2.50501e − 09 1.04889e − 09 6.77362e − 11 4.28567e − 12 1.10441e − 13
Table 3.2 From the above table, it shows that compact finite difference method can obtain error O(h5 ) at all the mesh points in the interval [a, b], but when h is too small, rounding errors take over. Last we will try different initial values vi0 , i = 1, · · · , N − 1, and see whether 1 , and we get the following table: our method depends on them or not. Let h = 50 tests test 1 test 2 test 3 test 4 test 5 test 6
Initial values vi0 vi0 = i/50, i = 1, · · · , 49 vi0 = 1 − i/50, i = 1, · · · , 49 vi0 = 0, i = 1, · · · , 49 vi0 = 100, i = 1, · · · , 49 vi0 = −8, i = 1, · · · , 49 vi0 = −10, i = 1, · · · , 49
Iteration numbers m 12 12 13 14 16 >100
Max absolute residual err. 1.50966e − 08 1.50966e − 08 1.50966e − 08 1.50966e − 08 1.50966e − 08 NaN
Table 3.3 It seems that the compact finite difference method converges very fast for test 1 ∼ test 5 and doesn’t depend too much on initial values. But for test 6, it doesn’t converge, the reason is that the term e−v(x) at the right side of the equation (1.2) is a large number when V 0 = −10. Suitably choosing V 0 can make us get desirable solutions V in a shorter time, i.e., smaller iterations m and guarantee this method converges. From the above numerical experiments, we can easily see that the compact finite difference method indeed obtain O(h5 ) and this method also converges fast. 3.5. Dealing with Different Boundary Conditions In this part, suppose we still have integro-differential equation (1.2), but with different boundary conditions instead of typical boundary conditions (1.3). When we have the following boundary conditions: u(a) = ua , u0 (b) = udb ,
(3.28)
u0 (a) = uda , u(b) = ub ,
(3.29)
or
we can deal with them almost in the same way as we did in section 2.2. The only change we need to make is that to use the following sixth order equation h2 127 00 86 00 257 00 461 00 0 (− vN −4 + vN vN −2 + v ) = vN −2 − vN −1 + hvN , (3.30) −3 − 12 30 5 10 15 N −1
14
Jichao Zhao and Robert M. Corless
for the differential boundary condition of (3.28) or h2 127 00 86 00 257 00 461 00 (− v + v3 − v + v ) = −v1 + v2 − hv00 , (3.31) 12 30 4 5 10 2 15 1 for the differential boundary condition of (3.29), instead of equation (2.6) and equation (2.5) respectively. Matrix M , A and vector H need to be modified accordingly and the algorithm is same with typical boundary conditions. In some conditions, we may meet with the following boundary conditions u(a) = ua , u0 (a) = uda ,
(3.32)
For this case, the boundary conditions are quite different with the boundary conditions we have mentioned before, since we only know the informations at one of two end points. The equation (2.6) can not be used since we don’t know the value of vN . We need one more equation except equation (2.5) and equations (2.1), so we try to use the known values, u0 (a) and u(a) together, to get another equation: h2 97 00 132 00 207 00 502 00 ( v − v + v − v ) = −2v1 + 2v0 + 2hv00 , (3.33) 12 15 4 5 3 5 2 15 1 also the above equation can obtain O(h6 ) order. Same thing for the following boundary conditions: u(b) = ub , u0 (b) = udb ,
(3.34)
We use the following equation except equation (2.1) and equation (2.6): h2 97 00 132 00 207 00 502 00 0 ( v − v + v − v ) = −2vN −1 + 2vN − 2hvN . (3.35) 12 15 N −4 5 N −3 5 N −2 15 N −1 For the following boundary conditions: u0 (a) = uda , u0 (b) = udb ,
(3.36)
we meet some problems. We may think of using equations similar to (3.31) and (3.30) to deal with boundary conditions (3.36), but no matter what accuracy of the formula will be, when we write them with equations (2.1), i = 1, · · · , N − 1 together into matrix form, compact finite difference matrix A always be: −1 1 0 ··· 0 0 .. .. 1 −2 1 . . . . . .. .. .. 0 0 0 A= , . .. .. .. .. . . . 0 0 1 −2 1 0 0 ··· 0 1 −1 but note det(A) = 0. So it means for boundary conditions (3.36), we need more conditions to solve integro-differential equation (1.2). In [1], it mentions the following general boundary conditions: a1 u(a) + b1 u(b) = λ1 , a2 u0 (a) + b2 u0 (b) = λ2 ,
(3.37)
Compact Finite Difference Method for IDE
15
where a1 , a2 , b1 , b2 , λ1 , and λ2 are known constants. Using Taylor expansion, we can use a linear combination of u(a), u1 , u2 , u3 , u4 and u5 to replace u0 (a) in (3.37) with sixth order of accuracy. Same way to replace u0 (b) by a linear combination of u(b), uN −1 , uN −2 , uN −3 , uN −4 and uN −5 . Solving for u(a) and u(b), we can get: u(a) = φ1 (u1 , u2 , u3 , u4 , u5 , uN −5 , uN −4 , uN −3 , uN −2 , uN −1 ), u(b) = φ2 (u1 , u2 , u3 , u4 , u4 , uN −5 , uN −4 , uN −3 , uN −2 , uN −1 ). Since we have the above relations, we can just deal with them in the same way as the typical boundary conditions by compact finite difference method. At the end of this part, we will give a procedure to show how to handle the case when we have the conditions are not limited on boundary points. Take a simple example to explain, suppose we have the following conditions: u(a) = ua , u(c) = uc ,
(3.38)
where ua and uc are two already known constants, and c ∈ (a, b). The procedure to handle them are as follows: 1. Solve the integro-differential equation on the interval [a, c] by compact finite difference method since we know the boundary conditions (3.38). So we can obtain all the numerical values for the mesh points in the subinterval [a, c]. 2. Construct another boundary condition u0 (c) by the linear combination of five mesh points in interval [a, c] closer to point c with the accuracy of sixth order by Taylor expansion. 3. Use the values of u0 (c) and u(c) as the boundary conditions to solve the integro-differential equation on the interval [c, b] by compact finite difference method. Then we obtain all the numerical values for the whole interval [a, b] with the condition (3.38). We can solve for other types of mixed boundary and interior points in the similar way.
4. Solving First Order Integro-Differential Equations We consider how to use compact finite difference method to solve the following first order integro−differential equation: Z b du(x) + K(x, s)u(s)ds = f (x, u), a ≤ x ≤ b, (4.1) − dx a with boundary condition u(a) = u0 .
(4.2)
The standard compact finite difference formula for first derivatives is: 0 0 h(vi−1 + 4vi0 + vi+1 ) = 3(−vi−1 + vi+1 ),
(4.3)
16
Jichao Zhao and Robert M. Corless
and the error in this formula is O(h5 ). To consistent with the truncation error of O(h6 ) for second derivatives in section 3 and get a better property of matrix A, we use the following formula for first derivatives: h (−23v10 + 39v20 − 21v30 + 5v40 ) = v0 − 2v1 + v2 , (4.4) 12 9 1 h (− v10 − 5v20 + 12v30 − 3v40 + v50 ) = v1 − 2v2 + v2 , (4.5) 12 2 2 h 1 0 1 0 0 0 ( v − 7vi−1 + 7vi+1 − vi+2 ) = vi−1 − 2vi + vi+1 , i = 3, · · · , N − 2, (4.6) 12 2 i−2 2 h 1 0 9 0 0 0 0 (− vN −4 + 3vN −3 − 12vN −2 + 5vN −1 + vN ) = vN −2 − 2vN −1 + vN , (4.7) 12 2 2 h 1 0 5 0 19 0 9 0 (− v + v − vN −1 − vN ) = vN −1 − vN , (4.8) 12 2 N −3 2 N −2 2 2 6 all formulae are O(h ), and can be obtained by Taylor expansion. Or rewrite them in matrix form: M V 0 = AV + H, (4.9) where A=
−2
1
1
−2 .. .
0 .. . 0
···
−23 9 − 2 0 . h .. M= 12 0 . .. ··· 0
H
0 0 v1 .. v20 1 . .. .. and V 0 = ... , . . 0 v 0 N −1 1 −2 1 0 vN 0 1 −1 N ×N 39 −21 5 0 ··· 0 .. −5 12 −3 − 21 0 . 1 1 −7 0 7 −2 0 2 .. .. .. .. .. .. . . . . . . , 1 1 −7 0 7 −2 0 2 1 0 −7 0 7 − 21 2 9 0 − 12 3 −12 5 2 ··· 0 −5 21 −39 23 N ×N v1 v0 v2 0 .. = . and V = ... . vN −1 0 vN 0 0
···
Compact Finite Difference Method for IDE Matrix A can 1 −1 2 L= 0 .. . 0
be factored exactly into LU : 0 0 ··· 0 −2 1 .. −3 0 1 0 . 2 .. .. .. . .. . . . 0 , U = 0 . −(N −2) .. 1 0 N −1 −(N −1) 0 ··· ··· 0 1 N
17
0
···
1 .. .
..
0
−N N −1
.
0
0 .. . . 0 1 −1 N
For the integral part of (4.1), we use the same way as we did for second order integro-differential equation. Then we can solve the first order integro−differential problem (4.1) and (4.2) with high order of accuracy by compact finite difference method. Finally we give an example to show our compact finite difference method works for first order integro−differential equation. Example 3: Z 1 du(x) − + e(−s) sin xu(s)ds = u(x) − 2x, 0 ≤ x ≤ 1, (4.10) dx 0 u(0) = 0, The kernel function K(x, s) = sin(x)e(−s) , f (x, u) = u(x) − 2x, and the exact solution is: u(x) = C1 sin(x) + 2x − 2 + C2 cos(x) + C3 e(−x) , C1 = .0686 · · · , C2 = −.0686 · · · , C3 = 2.0686 · · · , By compact finite difference method, we get the Fig. 3 and 4 for example 3.
5. Solving the System of Integro-Differential Equations When we have a system of integro-differential equations, we can also use compact finite differene method to solve them. To make life easier, we take J, the number of integro-differential equations, as the two and both of them are second order integro-differential equations, to explain how to solve the system. In the similar way, we can solve for the case J > 2 or the system has first order integro-differential equations. In this section, we will focus on the following system of two second order integro-differential equations: Z b Z b d2 u1 (x) − + K11 (x, s)u1 (s)ds + K12 (x, s)u2 (s)ds = f1 (x, u1 , u2 ), (5.1) dx a a Z Z b b d2 u2 (x) + K21 (x, s)u1 (s)ds + K22 (x, s)u2 (s)ds = f2 (x, u1 , u2 ), (5.2) − dx a a a ≤ x ≤ b,
18
Jichao Zhao and Robert M. Corless
7e-07
4e-07
Error 2e-07
1e-07 0.2
0.4
0.6
0.8
x Exact error Residual error
Figure 3. Exact and residual err. for example 3, h = .1 with boundary conditions u1 (a) = u1a , u1 (b) = u1b , u2 (a) = u2a , u2 (a) = u2b ,
(5.3) (5.4)
And the algorithm to solve the above (5.1∼5.4) is given below: 1. By compact finite difference method, we can write the differential parts of equations (5.1) and (5.2) into matrix form: M V100 = AV1 + H1 , M V200 = AV2 + H2 , where the elements of vector V1 and V2 correspond to partition points of functions u1 (x) and u2 (x), please refer to equation (2.7) of section 2 for details. 2. We can obtain the two following matrix equations after we combine the above two matrix equations and Gauss 3 point rule for the integral parts of equations (5.1) and (5.2): AV1 = G(V1 , V2 ), AV1 = G(V1 , V2 ).
(5.5)
Compact Finite Difference Method for IDE
19
1e-06
1e-07
1e-08
Error
1e-09
1e-10
1e-11
.5e-2
.1e-1
.2e-1
.5e-1
.1
h Exact error Residual error
Figure 4. Exact and residual err. with varying h for example 3
3. Set initial values for vector V10 and V20 , then compute new V1m+1 and V2m+1 using already known values of V1m and V2m by equations (5.5), m = 0, 1, · · · , until the difference, both V1m+1 − V1m and V2m+1 − V2m , is smaller than a given error. The following example, a system of integro-differential equations, is solved by compact finite difference method. Example 4: d2 u1 (x) − + dx2 −
Z
Z
1
sin (x)su1 (s)ds + 2 0
d2 u2 (x) + dx2
Z
Z
1
su1 (s)ds + 0
1
e−s u2 (s)ds = −u1 (x) − e−x ,
0 1
cos (x)e−s u2 (s)ds = −u2 (x) − 4,
0
0 ≤ x ≤ 1,
20
Jichao Zhao and Robert M. Corless
with boundary conditions u1 (0) = 0, u1 (1) = 1, u2 (0) = 0, u2 (1) = 1,
(5.6) (5.7)
and the analytic solution is: 1 u1 (x) = C10 + C11 sin(x) + C12 e(x) + C13 e(−x) − xe(−x) , 2 u2 (x) = C20 + C21 cos(x) + C22 e(x) + C23 e(−x) , C10 = −.0763 · · · , C11 = −.1410 · · · , C12 = .5747 · · · , C13 = −.4983 · · · , C20 = −4.2820 · · · , C21 = −.0190 · · · , C22 = 1.5784 · · · , C23 = 2.7226 · · · . By compact finite difference method, we get the numerical solutions for u1 and u2 , please see Fig. 5 ∼ 8 for them.
6. Conclusions Compact finite difference method can solve very general integro-differential equations with high accuracy, while only use time complexity O(N ) for solving the matrix equations. It is easy to extend our work for Volterra integro-differential equation.
References [1] A. Karameta, and M. Sezer, A Taylor collocation method for the solution of linear integro-differential equations, Intern. J. Computer Math., 2002, Vol. 79(9), pp. 9871000. [2] E. Liz, and J.J. Nieto, Boundary value problems for second order integro−differential equations of Fredholm type, J. Comput. Appl. Math. 72 (1996) 215−225. [3] J. Zhao, and T. Zhang, A highly accurate derivative recovery formula to integrodifferential equations, Numerical Mathematics a Journal of Chinese Universities, 2004 Vol.26(1), pp.81-90. [4] H. Sun, and J. Zhang, A high order compact boundary value method for solving one dimensional heat equations, Technical Report No. 333-02, University of Kentucky, 2002. [5] H.Yang, Numerical analysis of American options, Ph.D. Thesis, University of Alberta, 2001. [6] K. Sty´s, and T. Sty´s, A higher-order finite difference method for solving a system of integro-differential equations, Journal of Computational and Applied Mathematics, 126(2000) 33-46. [7] M.O. Ahmed, An exploration of compact finite difference methods for the numerical solution of PDE, Ph.D. thesis, the University of Western Ontario, 1997.
Compact Finite Difference Method for IDE
21
1e-06
8e-07
6e-07 Error
4e-07
3e-07 0.2
0.4
0.6
0.8
x Exact error Residual error
Figure 5. Exact and residual err. for u1 (x), h = 0.1
[8] R. Cont, P. Tankov and E. Voltchkova, Option pricing models with jumps: integrodifferential equations and inverse problems. European Congress on Computational Methods in Applied Sciences and Engineering, 24-28 July (2004) 1-20. [9] R.H. Li, Z.Y. Chen and W. Wu, Generalized difference methods for differential equations (Numerical analysis of finite volume methods), Marcel Dekker, Inc., New York, 2000. [10] R.S. Varga, Matrix iterative analysis, Berlin; New York: Springer, c2000. [11] R.M. Corless, and J. Rokicki, The symbolic generation of finite-difference formulas. Zeitschrift fur Angewandte Mathematik und Mechanik (ZAMM) Proceedings. ICIAM/GAMM 95, Numerical Analysis, Scientific Computing, Computer Science. Hamburg, Germany: 3 July 1995-7 July 1995. G. Alefeld, O. Mahrenholtz, and R. Mennicken. ZAMM, Vol. 76, Supplement. Germany: Akademie Verlag. 381-382.
22
Jichao Zhao and Robert M. Corless
2e-06
1e-06
Error 7e-07
4e-07
0.2
0.4
0.6
0.8
x Exact error Residual error
Figure 6. Exact and residual err. for u2 (x), h = 0.1
[12] R.M. Corless, J. Rokicki and J. Zhao, FINDIF, a routine for generation of finite difference formulae, Share Library package 1994, and upgraded to n dimensions for “iguana”. [13] R. Natalini, M. Briani and C.L. Chioma, Finite differences schemes for integrodifferential equations in financial mathematics, IAC-CNR & University, 2003. [14] T. Zhang, Finite element methods for evolution integro-differential equations, Northeastern University Press, Shenyang of China, 2001. [15] W. Liu, The exponential stabilization of the high dimensional linear system of thermoviscoelasticity. J. Math. Pures Appl., 77 (1998) 355-386.
Compact Finite Difference Method for IDE
23
1e-06
1e-07
1e-08 Error
1e-09
1e-10
1e-11 .5e-2
.1e-1
.2e-1 h
.5e-1 Exact error Residual error
Figure 7. Exact and residual err. for u1 (x) with varying h Acknowledgment Zhao is thankful to Prof. T. Zhang for remarks on a draft of this paper. Jichao Zhao Department of Applied Mathematics, University of Western Ontario, Middlesex College, London, Ontario, Canada N6A 5B7 e-mail: [email protected] Robert M. Corless Department of Applied Mathematics, University of Western Ontario, Middlesex College, London, Ontario, Canada N6A 5B7 e-mail: [email protected]
.1
24
Jichao Zhao and Robert M. Corless
1e-06
1e-07 1e-08 Error 1e-09
1e-10 1e-11
.5e-2
.1e-1
.2e-1 h
.5e-1 Exact error Residual error
Figure 8. Exact and residual err. for u2 (x) with varying h
.1
