View source for Orthogonal double-sweep method

<!--
o0702901.png
$#A+1 = 53 n = 0
$#C+1 = 53 : ~/encyclopedia/old_files/data/O070/O.0700290 Orthogonal double\AAhsweep method
Automatically converted into TeX, above some diagnostics.
Please remove this comment and the {{TEX|auto}} line below,
if TeX found to be correct.
-->

{{TEX|auto}}
{{TEX|done}}

A variant of the [[Double-sweep method|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 [[#References|[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 [[#References|[2]]]).

Methods based on the use of a [[Fundamental system of solutions|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 [[#References|[1]]], [[#References|[3]]]). However, these methods are really variants of the [[Shooting method|shooting method]].

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

====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|invariant imbedding]]), which however uses one non-linear equation only. Other variants, with similar ideas as this  "orthogonal double-sweep method" , are in [[#References|[a2]]].

====References====
<table><TR><TD valign="top">[a1]</TD> <TD valign="top">  U.M. Ascher,   R.M.M. Mattheij,   R.D. Russell,   "Numerical solution for boundary value problems for ordinary differential equations" , Prentice-Hall  (1988)</TD></TR><TR><TD valign="top">[a2]</TD> <TD valign="top">  G.M. Meyer,   "Continuous orthonormalization for boundary value problems"  ''J. Comput. Phys.'' , '''62'''  (1986)  pp. 248–262</TD></TR></table>