Orthogonal double-sweep method

A variant of the double-sweep method based on an orthogonal transformation of the unknowns. Let, for $ a \leq x \leq b $, a boundary value problem be examined for a pair of linear ordinary differential equations

$$ \tag{1 } y ^ \prime ( x) = a _ {1} ( x) y( x) + b _ {1} ( x) z( x) + f _ {1} ( x), $$

$$ \tag{2 } z ^ \prime ( x) = a _ {2} ( x) y( x) + b _ {2} ( x) z( x) + f _ {2} ( x), $$

with conditions of the form

$$ \tag{3 } \alpha _ {1} y( a) + \beta _ {1} z( a) = \gamma _ {1} ,\ \ \alpha _ {1} ^ {2} + \beta _ {1} ^ {2} = 1, $$

$$ \tag{4 } \alpha _ {2} y( b) + \beta _ {2} z( b) = \gamma _ {2} ,\ \ \alpha _ {2} ^ {2} + \beta _ {2} ^ {2} = 1. $$

Let the given functions $ a _ {i} ( x), b _ {i} ( x), f _ {i} ( x) $, $ i = 1, 2 $, be continuous on the segment $ a \leq x \leq b $. A solution of the boundary value problem (1)–(4) by the orthogonal double-sweep method is realized as follows.

I) The auxiliary Cauchy problem

$$ s ^ \prime ( x) = c( x) r( x),\ c ^ \prime ( x) = - s( x) r( x), $$

$$ r( x) = s ^ {2} ( x) b _ {1} ( x) + s( x) c( x)( b _ {2} ( x)- a _ {1} ( x)) + $$

$$ \tag{5 } - c ^ {2} ( x) a _ {2} ( x), $$

$$ \tag{6 } u ^ \prime ( x) = \overline{a}\; _ {1} ( x) u( x) + \overline{f}\; _ {1} ( x), $$

$$ \tag{7 } s( a) = \alpha _ {1} ,\ c( a) = \beta _ {1} ,\ u( a) = \gamma _ {1} , $$

is solved, where

$$ \overline{a}\; _ {1} ( x) = a _ {1} ( x) s ^ {2} ( x) + b _ {2} ( x) c ^ {2} ( x) + $$

$$ + ( b _ {1} ( x) + a _ {2} ( x)) s( x) c( x), $$

$$ \overline{f}\; _ {1} ( x) = f _ {1} ( x) s( x) + f _ {2} ( x) c( x) $$

(the direct double-sweep).

II) The condition $ \Delta = \alpha _ {2} c( b) - \beta _ {2} s( b) \neq 0 $ is tested, and if it is fulfilled, the Cauchy problem

$$ \tag{8 } v ^ \prime ( x) = \overline{a}\; _ {2} ( x) u( x) + \overline{b}\; _ {2} ( x) v( x) + \overline{f}\; _ {2} ( x), $$

$$ \tag{9 } v( b) = \frac{\gamma _ {2} - [ \alpha _ {2} s( b) + \beta _ {2} c( b)] u( b) } \Delta , $$

where

$$ \overline{a}\; _ {2} ( x) = 2( a _ {1} ( x) - b _ {2} ( x)) s( x) c( x) + $$

$$ + ( b _ {1} ( x) + a _ {2} ( x))( c ^ {2} ( x) - s ^ {2} ( x)), $$

$$ \overline{b}\; _ {2} ( x) = a _ {1} ( x) c ^ {2} ( x) + b _ {2} ( x) s ^ {2} ( x) + $$

$$ - ( b _ {1} ( x) + a _ {2} ( x)) s( x) c( x), $$

$$ \overline{f}\; _ {2} ( x) = f _ {1} ( x) c( x) - f _ {2} ( x) s( x), $$

is solved in the direction from $ x = a $ to $ x = b $( the inverse double-sweep).

III) The required functions are calculated using the formulas

$$ y( x) = s( x) u( x) + c( x) v( x), $$

$$ z( x) = c( x) u( x) - s( x) v( x). $$

If the solution $ y( x), z( x) $ of the boundary value problem (1)–(4) exists and is unique and stable with respect to small changes of the coefficients and the free terms defining the problem, then $ \Delta \neq 0 $ and the method in question is also stable (see [2]).

A system of linear algebraic equations

$$ \tag{10 } y _ {k+} 1 = A _ {k} y _ {k} + B _ {k} z _ {k} + F _ {k} , $$

$$ \tag{11 } z _ {k+} 1 = C _ {k} y _ {k} + D _ {k} z _ {k} + G _ {k} ,\ k = 0 \dots n- 1, $$

$$ \tag{12 } \alpha _ {0} y _ {0} + \beta _ {0} z _ {0} = \gamma _ {0} , $$

$$ \tag{13 } \alpha _ {n} y _ {n} + \beta _ {n} z _ {n} = \gamma _ {n} , $$

where $ A _ {k} D _ {k} \neq B _ {k} C _ {k} $, $ \alpha _ {0} ^ {2} + \beta _ {0} ^ {2} = 1 $, $ \alpha _ {n} ^ {2} + \beta _ {n} ^ {2} = 1 $, is solved according to the following rules.

1) Using the formulas

$$ s _ {k+} 1 = \frac{C _ {k} c _ {k} - D _ {k} s _ {k} }{\rho _ {k} } , $$

$$ c _ {k+} 1 = \frac{B _ {k} s _ {k} - A _ {k} c _ {k} }{\rho _ {k} } , $$

$$ \rho _ {k} = \sqrt {( C _ {k} c _ {k} - D _ {k} s _ {k} ) ^ {2} + ( B _ {k} s _ {k} - A _ {k} c _ {k} ) ^ {2} } , $$

$$ u _ {k+} 1 = ( A _ {k} s _ {k} s _ {k+} 1 + B _ {k} c _ {k} s _ {k+} 1 + C _ {k} s _ {k} c _ {k+} 1 + D _ {k} c _ {k} c _ {k+} 1 ) u _ {k} + $$

$$ + ( F _ {k} s _ {k+} 1 + G _ {k} c _ {k+} 1 ), $$

$$ s _ {0} = \alpha _ {0} ,\ c _ {0} = \beta _ {0} ,\ u _ {0} = \gamma _ {0} , $$

the quantities $ s _ {k+} 1 , c _ {k+} 1 , u _ {k+} 1 $ are calculated successively when $ k = 0 \dots n- 1 $( the direct double-sweep).

2) The condition $ \Delta _ {n} = \alpha _ {n} c _ {n} - \beta _ {n} s _ {n} \neq 0 $ is tested, and if it is fulfilled, then

$$ v _ {n} = \frac{\gamma _ {n} - ( \alpha _ {n} s _ {n} + \beta _ {n} c _ {n} ) u _ {n} }{\Delta _ {n} } $$

and

$$ v _ {k} = \frac{1}{\rho _ {k} } \{ v _ {k+} 1 + $$

$$ + [( C _ {k} s _ {k} + D _ {k} c _ {k} ) s _ {k+} 1 - ( A _ {k} s _ {k} + B _ {k} c _ {k} ) c _ {k+} 1 ] u _ {k} + $$

$$ + {} ( G _ {k} s _ {k+} 1 - F _ {k} c _ {k+} 1 ) \} $$

are calculated, when $ k = n- 1 , n- 2 \dots 1 $( inverse double-sweep).

3) The values of the required solution of the system of equations (10)–(13) are calculated using the formulas

$$ y _ {k} = u _ {k} s _ {k} + v _ {k} c _ {k} ,\ \ z _ {k} = u _ {k} c _ {k} - v _ {k} s _ {k} . $$

If a solution of the system of equations (10)–(13) exists and is unique and stable with respect to small changes of the coefficients and the free terms, then the orthogonal double-sweep method in question is also stable (see [2]).

Methods based on the use of a fundamental system of solutions of a homogeneous system of equations which aim at transferring the boundary conditions are sometimes called orthogonal double-sweep methods (see [1], [3]). However, these methods are really variants of the shooting method.

References

[1]	N.S. Bakhvalov, "Numerical methods: analysis, algebra, ordinary differential equations" , MIR (1977) (Translated from Russian)
[2]	V.I. Krylov, V.V. Bobkov, P.I. Monastyrnyi, "Computing methods of higher mathematics" , 2 , Minsk (1975) (In Russian)
[3]	A.A. Samarskii, E.S. Nikolaev, "Numerical methods for grid equations" , 1–2 , Birkhäuser (1989) (Translated from Russian)

Comments

Care should be taken with some steps in the above algorithm, as severe cancellation may occur (e.g. in $ \Delta $). Also, the method requires solving non-linear equations, whereas the system is simple and linear. This aspect is also shared by the Riccati method (or invariant imbedding), which however uses one non-linear equation only. Other variants, with similar ideas as this "orthogonal double-sweep method" , are in [a2].

References

[a1]	U.M. Ascher, R.M.M. Mattheij, R.D. Russell, "Numerical solution for boundary value problems for ordinary differential equations" , Prentice-Hall (1988)
[a2]	G.M. Meyer, "Continuous orthonormalization for boundary value problems" J. Comput. Phys. , 62 (1986) pp. 248–262