Accurate numerical solutions to stationary free surface problems from capillarity

Abstract

The paper considers two static problems from capillarity. The first one consists in the determination of the surface of a liquid in a capillary tube and the second in the computation of the shape of a sessile or pendent drop of liquid of a given volume, both configurations being considered in the gravitational field. From the mathematical point of view both problems are nonlinear two-point boundary value problems which include in their formulations additional geometrical unknowns, i.e., the so-called free boundary value problems. Computationally they are laborious problems since they involve non-normal Jacobian matrices. The resulting numerical difficulties are considerably reduced by treating these problems by a collocation method in conjunction with the continuation with respect to the parameters. The method is implemented by the MATLAB code bvp4c. The numerical evaluations of the free surfaces are in reasonable accordance with asymptotic estimations.

Authors

Călin-Ioan Gheorghiu
(Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy)

Keywords

capillarity; free boundary value problem; liquid gas interface; sessile; pendent drop; collocation; matlab; bvp4c.

References

See the expanding block below.

Cite this paper as

C.I. Gheorghiu, Accurate numerical solutions to stationary free surface problems from capillarity, Phys. Part. Nuclei Lett., 5 (2008), 286-289.
doi: 10.1134/s154747710803031x

PDF

About this paper

Journal

Physics of Particles and Nuclei Letters

Publisher Name

SP MAIK Nauka/Interperiodica

Print ISSN

1547-4771

Online ISSN

1531-8567

Google Scholar Profile

?

[1] Finn R. Capillary Free Surfaces. Preprint No. 7/1981. Bonn: Inst. féur Angewante Mathematik Univ.,1981.
Google Scholar

2. Gheorghiu C.-I., Mathematica. Studia Univ. Babes-Bolyai. 2005. V. L. P. 61.

3. Kierzenka J., Shampine L. F. // ACM Trans. Math. Softw. 2001. V. 27. P. 299.

4. Miele A., Aggarwal A. K., Tietze J. L. // J. Comp. Phys. 1974. V. 15. P. 117.

5. Trefethen L. N., Bau III D. Numerical Linear Algebra. SIAM, 1997.

Paper (preprint) in HTML form

S154747710803031X

Accurate Numerical Solutions to Stationary Free Surface Problems from Capillarity 1 1 ^(1){ }^{\mathbf{1}}1

C.-I. Gheorghiu"T. Popoviciu" Institute of Numerical Analysis, Cluj-Napoca, Romania

Abstract

The paper considers two static problems from capillarity. The first one consists in the determination of the surface of a liquid in a capillary tube and the second in the computation of the shape of a sessile or pendent drop of liquid of a given volume, both configurations being considered in the gravitational field. From a mathematical point of view, both problems are nonlinear two-point boundary value problems which include in their formulations additional geometrical unknowns, i.e., the so-called free boundary value problems. Computational, they are laborious problems owing to the fact that they involve nonnormal Jacobian matrices. The resulting numerical difficulties are considerably reduced by treating these problems by a collocation method in conjunction with the continuation with respect to the parameters. The method is implemented by the MATLAB code b v p 4 c b v p 4 c bvp4cb v p 4 cbvp4c. The numerical evaluations of the free surfaces are in reasonable accordance with asymptotic estimations.

INTRODUCTION

In many cases, problems in viscous fluid mechanics can be described by boundary value problems for partial differential equations. In some cases, using the geometrical symmetries of the configurations, such problems can be reduced to two-point boundary value problems. When there are additional geometrical unknowns, we speak of free boundary value problems (FB for short). More specific, these geometrical unknowns can be some values of the unknowns on subsets of boundary points or even the coordinates of the boundary points themselves. In this paper we solve numerically two classes of static capillary FBs. These problems are defined as FBs on which we must take into account capillary forces, i.e., forces due to the intermolecular attractions, which have a nonvanishing resultant at the boundary of the liquid.
The first one is to determine the liquid-air interface (surface) whose mean curvature is proportional to its height above the reference plane, and which makes with a container bounding wall a prescribed constant angle. The second one is to find the shape of a sessile or pendent liquid drop of a given volume, which rests in the gravitational field. The actual computations were carried out for symmetric cases. The equations of these problems, as well as the analytic properties of the solutions (existence, uniqueness, and asymptotic behavior for Bond number approaching 0), were established in the monograph of Finn [1]. Both problems lead to nonlinear second order two-point free boundary value problems, i.e., the values of the solution at some bound-
ary points are unknown, but some values of the derivatives are supplied in order to close the problem. Their solutions depend on the nondimensional Bond number (Bo for short) which is a measure of the "size" of the configuration, i.e., it is the ratio of gravitational forces to the capillary (superficial) forces. The main computational difficulty with these problems consists in the fact that they involve nonnormal Jacobian matrices. A test problem with a Jacobian matrix in close form is considered in order to underline these aspects.
The problems are solved numerically making use of a collocation type method, based on cubic interpolation and introduced in the paper of Kierzenka and Shampine [3]. It is implemented by the MATLAB code b v p 4 c b v p 4 c bvp4cb v p 4 cbvp4c and provides a smooth solution continuous along with its first derivative. In contrast with well-known shooting, this solution is approximated on the whole interval and the boundary conditions are taken into account at all times. As the resulting nonlinear system of algebraic equations is solved iteratively, the question of coming up with a sufficiently accurate initial guess of unknowns is particularly discussed.

1. A TEST PROBLEM

In their paper [4], Miele, Aggarwal, and Tietze consider the following two-point boundary value problem introduced by Troesch, namely,
(1) { x t = y y t = k sinh ( k x ) , 0 < t < 1 x ( 0 ) = 0 , x ( 1 ) = 1 . (1) x t = y y t = k sinh ( k x ) , 0 < t < 1 x ( 0 ) = 0 , x ( 1 ) = 1 . {:(1){[x_(t)=y],[y_(t)=k sinh(kx)","quad0 < t < 1],[x(0)=0","quad x(1)=1.]:}:}\left\{\begin{array}{l} x_{t}=y \tag{1}\\ y_{t}=k \sinh (k x), \quad 0<t<1 \\ x(0)=0, \quad x(1)=1 . \end{array}\right.(1){xt=yyt=ksinh(kx),0<t<1x(0)=0,x(1)=1.
The Jacobian matrix,
J = ( 0 k 2 cosh ( k x ) 1 0 ) J = 0 k 2 cosh ( k x ) 1 0 J=([0,k^(2)cosh(kx)],[1,0])J=\left(\begin{array}{cc} 0 & k^{2} \cosh (k x) \\ 1 & 0 \end{array}\right)J=(0k2cosh(kx)10)
is characterized by the eigenvalues λ = ± k cosh ( k x ) λ = ± k cosh ( k x ) lambda=+-ksqrt(cosh(kx))\lambda= \pm k \sqrt{\cosh (k x)}λ=±kcosh(kx), which at the endpoints become λ ( 0 ) = ± k λ ( 0 ) = ± k lambda(0)=+-k\lambda(0)= \pm kλ(0)=±k and λ ( 1 ) = ± k cosh ( k ) λ ( 1 ) = ± k cosh ( k ) lambda(1)=+-ksqrt(cosh(k))\lambda(1)= \pm k \sqrt{\cosh (k)}λ(1)=±kcosh(k).
For relatively low values of k k kkk, i.e., 0 < k 4 0 < k 4 0 < k <= 40<k \leq 40<k4, the above eigenvalues remain relatively small and the problem can be solved employing the usual finite difference or the shooting methods. On the other hand, for relatively large values of k k kkk, i.e., 5 k 10 5 k 10 5 <= k <= 105 \leq k \leq 105k10, the above quoted authors make use of a special multipoint boundaryvalue technique in order to solve the problem. They observe that the problem is characterized by the first integral
(2) z := cosh ( k x ) y 2 / 2 = const. (2) z := cosh ( k x ) y 2 / 2 =  const.  {:(2)z:=cosh(kx)-y^(2)//2=" const. ":}\begin{equation*} z:=\cosh (k x)-y^{2} / 2=\text { const. } \tag{2} \end{equation*}(2)z:=cosh(kx)y2/2= const. 
and use that as a "tool" in order to check the accuracy of the numerical solution. In our numerical experiments performed on a Pentium IV, 3.0 G , PC, the constancy of z z zzz is verified to a maximal error equal to 4.082997 e 06 4.082997 e 06 4.082997 e-064.082997 e-064.082997e06 for k = 5 k = 5 k=5k=5k=5, to an error which equals 7.651421e -03 for k = 10 k = 10 k=10k=10k=10, and to an error equal to 6.602193 e + 00 6.602193 e + 00 6.602193 e+006.602193 e+006.602193e+00 corresponding to k = 15 k = 15 k=15k=15k=15. Our numerical results are with one significant figure worse than those reported in [4, Table 2, p. 126], for k = 5 k = 5 k=5k=5k=5 and 10, respectively. For k k kkk larger than 16 , the accuracy in the conservation of z z zzz is lost rapidly. In spite of this fact, even at these values of k k kkk, our numerical solutions behave correctly in the left neighborhood of t = 1 t = 1 t=1t=1t=1. They show a boundary layer type behavior in this region.
It is worth noting at this point that we used essentially the technique of continuation with respect to parameter k k kkk in order to handle the problem for k 15 k 15 k >= 15k \geq 15k15. Without this technique for such values of k k kkk, the collocation equations become unsolvable due to singular Jacobian.
Remark 1. In fact the above Jacobian matrix J J JJJ is nonnormal, i.e., J J J J 0 J J J J 0 J^(**)J-JJ^(**)!=0J^{*} J-J J^{*} \neq 0JJJJ0, where *** stands for the conjugate transpose of a matrix (see Trefethen [5, p. 187]). The pseudospectrum of the matrix J J JJJ, at t = 1 t = 1 t=1t=1t=1, i.e., the spectrum of the randomly perturbed matrix with an arbitrary small quantity, is depicted in Fig. 1. In our previous paper [2], we introduced a scalar measure of nonnormality. It reads H ( J ) := ε ( J J J J ) / ε ( J ) H ( J ) := ε J J J J / ε ( J ) H(J):=sqrt(epsi(J^(**)J-JJ^(**)))//epsi(J)H(J):=\sqrt{\varepsilon\left(J^{*} J-J J^{*}\right)} / \varepsilon(J)H(J):=ε(JJJJ)/ε(J) and we noticed that 0 H ( J ) 2 4 , 0 0 H ( J ) 2 4 , 0 0 <= H(J) <= root(4)(2),00 \leq H(J) \leq \sqrt[4]{2}, 00H(J)24,0 being attained just in case of normal (symmetric) matrices. On this scale our matrix J J JJJ is situated extremely close to the upper bound. This scalar measure as well as the pseudospectrum lead us to the conclusion that the matrix J J JJJ is a genuine nonnormal one.
Fig. 1. The pseudospectrum of matrix J J JJJ for k = 4 k = 4 k=4k=4k=4.

2. THE LIQUID-GAS INTERFACE IN A CAPILLARY TUBE

In this section we consider static capillary liquidgas interfaces of liquids partly filling a vessel. More exactly we examine the behavior of capillary surfaces in a tube with circular section of radius R = a R = a R=aR=aR=a. We assume the following system of two coupled nonlinear differential equations and boundary conditions (see [1, p. 11])
(3) { div ( Tu ) = Bo u , in Ω v Tu = β , on Σ := Ω , (3) div ( Tu ) = Bo u ,  in  Ω v Tu = β ,  on  Σ := Ω , {:(3){[div(Tu)=Bou","quad" in "Omega],[vTu=beta","quad" on "quad Sigma:=del Omega","]:}:}\left\{\begin{array}{l} \operatorname{div}(\mathrm{Tu})=\mathrm{Bo} u, \quad \text { in } \Omega \tag{3}\\ v \mathrm{Tu}=\beta, \quad \text { on } \quad \Sigma:=\partial \Omega, \end{array}\right.(3){div(Tu)=Bou, in ΩvTu=β, on Σ:=Ω,
where, u ( x , y ) u ( x , y ) u(x,y)u(x, y)u(x,y) denotes the height over the wetted section Ω R 2 , T ( u ) := u 1 + | u | 2 , β Ω R 2 , T ( u ) := u 1 + | u | 2 , β Omega subR^(2),T(u):=(grad u)/(sqrt(1+|grad u|^(2))),beta\Omega \subset \mathbb{R}^{2}, T(u):=\frac{\nabla u}{\sqrt{1+|\nabla u|^{2}}}, \betaΩR2,T(u):=u1+|u|2,β is known, and v v vvv is the outward normal to Σ Σ Sigma\SigmaΣ.
Due to the symmetry with respect to the axis of the tube, on the half of the axial section, i.e., 0 r 1 , r := R / a 0 r 1 , r := R / a 0 <= r <= 1,r:=R//a0 \leq r \leq 1, r:= R / a0r1,r:=R/a, the problem reads
(4) { ( r u r 1 + u r 2 ) r = Bo r u , 0 < r < 1 u ( 0 ) = p 1 , u r ( 0 ) = 0 u ( 1 ) = p 2 , cos ( δ ) = β , r = 1 (4) r u r 1 + u r 2 r =  Bo  r u , 0 < r < 1 u ( 0 ) = p 1 , u r ( 0 ) = 0 u ( 1 ) = p 2 , cos ( δ ) = β , r = 1 {:(4){[((ru_(r))/(sqrt(1+u_(r)^(2))))_(r)=" Bo "ru","quad0 < r < 1],[u(0)=p_(1)","quadu_(r)(0)=0],[u(1)=p_(2)","quad cos(delta)=beta","quad r=1]:}:}\left\{\begin{array}{l} \left(\frac{r u_{r}}{\sqrt{1+u_{r}^{2}}}\right)_{r}=\text { Bo } r u, \quad 0<r<1 \tag{4}\\ u(0)=p_{1}, \quad u_{r}(0)=0 \\ u(1)=p_{2}, \quad \cos (\delta)=\beta, \quad r=1 \end{array}\right.(4){(rur1+ur2)r= Bo ru,0<r<1u(0)=p1,ur(0)=0u(1)=p2,cos(δ)=β,r=1
where Bo := ρ g a 2 / σ , ρ := ρ g a 2 / σ , ρ :=rho ga^(2)//sigma,rho:=\rho g a^{2} / \sigma, \rho:=ρga2/σ,ρ and σ σ sigma\sigmaσ being respectively the density and the surface tension of the liquid, and g g ggg is the usual gravitational acceleration. The parameters p 1 p 1 p_(1)p_{1}p1 and p 2 p 2 p_(2)p_{2}p2 are unknown and signify the heights of the liquid on the axis of symmetry and on the wall of the tube, respectively. The boundary condition u r ( 0 ) = 0 u r ( 0 ) = 0 u_(r)(0)=0u_{r}(0)=0ur(0)=0 appears
Fig. 2. The liquid-gas interface.
due to the symmetry and δ δ delta\deltaδ is the angle between the liq-uid-gas interface u ( r ) u ( r ) u(r)u(r)u(r) and the wall of the tube (see Fig. 2). The problem is effectively integrated as a first order differential system, i.e.,
(5) { u r = v v r = 1 r ( v + v 3 ) + Bo u ( v + v 3 ) 3 / 2 , 0 < r < 1 u ( 0 ) = p 1 , u ( 1 ) = p 2 v ( 0 ) = 0 , v ( 1 ) 1 + v 2 ( 1 ) = β . (5) u r = v v r = 1 r v + v 3 + Bo u v + v 3 3 / 2 , 0 < r < 1 u ( 0 ) = p 1 , u ( 1 ) = p 2 v ( 0 ) = 0 , v ( 1 ) 1 + v 2 ( 1 ) = β . {:(5){[u_(r)=v],[v_(r)=-(1)/(r)(v+v^(3))+Bo u(v+v^(3))^(3//2)","quad0 < r < 1],[u(0)=p_(1)","quad u(1)=p_(2)],[v(0)=0","quad(v(1))/(sqrt(1+v^(2)(1)))=beta.]:}:}\left\{\begin{array}{l} u_{r}=v \tag{5}\\ v_{r}=-\frac{1}{r}\left(v+v^{3}\right)+\operatorname{Bo} u\left(v+v^{3}\right)^{3 / 2}, \quad 0<r<1 \\ u(0)=p_{1}, \quad u(1)=p_{2} \\ v(0)=0, \quad \frac{v(1)}{\sqrt{1+v^{2}(1)}}=\beta . \end{array}\right.(5){ur=vvr=1r(v+v3)+Bou(v+v3)3/2,0<r<1u(0)=p1,u(1)=p2v(0)=0,v(1)1+v2(1)=β.
The second differential equation in (5) is singular at r = 0 r = 0 r=0r=0r=0, but the singular coefficient arises from the coordinate system and we expect a smooth solution in the neighborhood of this point. We deal with this singularity using the technique suggested in Example 6 from [3, pp. 18-19]. It means that we let r 0 r 0 r longrightarrow0r \longrightarrow 0r0 in this equation and obtain the value that must be used for evaluating the differential equation when r = 0 r = 0 r=0r=0r=0.
The existence of symmetric solutions with prescribed angle δ δ delta\deltaδ has been proved and the uniqueness theorem shows that the symmetric solutions are the only ones (see [1] and the paper quoted there). For Bo 0 Bo 0 Bolongrightarrow0\mathrm{Bo} \longrightarrow 0Bo0, in the same monograph, the Laplace formula, which estimates p 1 p 1 p_(1)p_{1}p1, and estimations for the parameter p 2 p 2 p_(2)p_{2}p2 are also available. These estimations were fairly useful in guessing the initial shape of the unknown u ( r ) u ( r ) u(r)u(r)u(r).
For Bo = 0.075 Bo = 0.075 Bo=0.075\mathrm{Bo}=0.075Bo=0.075 and Bo = 10 Bo = 10 Bo=10\mathrm{Bo}=10Bo=10, the dependence of u ( r ) u ( r ) u(r)u(r)u(r) on the angle δ δ delta\deltaδ is depicted in Fig. 2.
Remark 2. From the computational point of view, it is suggestive to notice that the code bvp4c works such that the maximal residual belongs to the range 4.563 e 4.563 e 4.563 e-4.563 e-4.563e 007 and 7.354 e 004 7.354 e 004 7.354 e-0047.354 e-0047.354e004.
Fig. 3. The shape of the sessile drops.

3. THE SESSILE AND PENDENT DROPS

We consider a liquid drop of volume V V VVV resting on a horizontal plane (surface) in a vertical gravity field and restrict our attention to the case of the contact angle δ δ delta\deltaδ such that 0 < δ π / 2 0 < δ π / 2 0 < delta <= pi//20<\delta \leq \pi / 20<δπ/2 (see Fig. 3). The case δ > π / 2 δ > π / 2 delta > pi//2\delta>\pi / 2δ>π/2 reduces to that one under the transformation u u u u u longrightarrow-uu \longrightarrow-uuu. Supposing that the drop is symmetric, the equations for the surface of the drop are
(6) { u r = v v r = 1 r ( v + v 3 ) + ( Bo u + λ ) ( v + v 3 ) 3 / 2 , 0 < r < 1 u ( 0 ) = p , v ( 0 ) = 0 u ( 1 ) = 0 , v ( 1 ) 1 + v 2 ( 1 ) = β , (6) u r = v v r = 1 r v + v 3 + ( Bo u + λ ) v + v 3 3 / 2 , 0 < r < 1 u ( 0 ) = p , v ( 0 ) = 0 u ( 1 ) = 0 , v ( 1 ) 1 + v 2 ( 1 ) = β , {:(6){[u_(r)=v],[v_(r)=-(1)/(r)(v+v^(3))+(Bou+lambda)(v+v^(3))^(3//2)","],[0 < r < 1],[u(0)=p","quad v(0)=0],[u(1)=0","quad(v(1))/(sqrt(1+v^(2)(1)))=beta","]:}:}\left\{\begin{array}{l} u_{r}=v \tag{6}\\ v_{r}=-\frac{1}{r}\left(v+v^{3}\right)+(\mathrm{Bo} u+\lambda)\left(v+v^{3}\right)^{3 / 2}, \\ 0<r<1 \\ u(0)=p, \quad v(0)=0 \\ u(1)=0, \quad \frac{v(1)}{\sqrt{1+v^{2}(1)}}=\beta, \end{array}\right.(6){ur=vvr=1r(v+v3)+(Bou+λ)(v+v3)3/2,0<r<1u(0)=p,v(0)=0u(1)=0,v(1)1+v2(1)=β,
where λ λ lambda\lambdaλ is a Lagrangian multiplier corresponding to the volume constraint. It is shown (see [1, pp. 59-60]) that, if δ = δ = delta=\delta=δ= const, then all equilibrium surfaces are symmetric and there is a one-to-one correspondence between capillary surfaces furnished by (4) and sessile drops implied by (6). Our numerical results are reported in Fig. 3.
Remark 3. An upper bound for the existence of the solution with respect to Bo can be obtained. The principle of virtual work applied to the energy functional of the drop of volume V V VVV yields a slightly modified form of (3), namely,
(7) { div ( Tu ) = Bo u + λ in Ω v Tu = β on Σ := Ω , (7) div ( Tu ) = Bo u + λ  in  Ω v Tu = β  on  Σ := Ω , {:(7){[div(Tu)=Bou+lambdaquad" in "Omega],[vTu=betaquad" on "Sigma:=del Omega","]:}:}\left\{\begin{array}{l} \operatorname{div}(\mathrm{Tu})=\mathrm{Bo} u+\lambda \quad \text { in } \Omega \tag{7}\\ v \mathrm{Tu}=\beta \quad \text { on } \Sigma:=\partial \Omega, \end{array}\right.(7){div(Tu)=Bou+λ in ΩvTu=β on Σ:=Ω,
where λ λ lambda\lambdaλ is associated with the volume restriction. The volume V V VVV of the drop can be expressed explicitly in terms of the boundary data and this property holds for general section Ω Ω Omega\OmegaΩ. We integrate the first equation in (3) and introduce the boundary condition. Thus, the volume is expressed as
(8) V = 1 Bo ( Σ cos δ λ Ω ) , (8) V = 1 Bo ( Σ cos δ λ Ω ) , {:(8)V=(1)/(Bo)(Sigma cos delta-lambda Omega)",":}\begin{equation*} V=\frac{1}{\mathrm{Bo}}(\Sigma \cos \delta-\lambda \Omega), \tag{8} \end{equation*}(8)V=1Bo(ΣcosδλΩ),
where the symbols Σ Σ Sigma\SigmaΣ and Ω Ω Omega\OmegaΩ denote, respectively, the measures of these sets. For symmetric solutions, (8) becomes
V = π Bo ( 2 cos δ λ ) V = π Bo ( 2 cos δ λ ) V=(pi)/(Bo)(2cos delta-lambda)V=\frac{\pi}{\mathrm{Bo}}(2 \cos \delta-\lambda)V=πBo(2cosδλ)
Taking into account the boundedness of cos δ cos δ cos delta\cos \deltacosδ, the last equality implies
(9) Bo 2 λ 2 0 1 u ( r ) d r (9) Bo 2 λ 2 0 1 u ( r ) d r {:(9)Bo <= (2-lambda)/(2int_(0)^(1)u(r)dr):}\begin{equation*} \mathrm{Bo} \leq \frac{2-\lambda}{2 \int_{0}^{1} u(r) d r} \tag{9} \end{equation*}(9)Bo2λ201u(r)dr
As the Lagrangian multiplier λ λ lambda\lambdaλ is delivered by b v p 4 c b v p 4 c bvp4cb v p 4 cbvp4c, this upper bound for Bo can be verified numerically. The results of our numerical experiments show that the quantity 2 λ 2 0 1 u ( r ) d r 2 λ 2 0 1 u ( r ) d r (2-lambda)/(2int_(0)^(1)u(r)dr)\frac{2-\lambda}{2 \int_{0}^{1} u(r) d r}2λ201u(r)dr overestimates the Bo number for
Bo 10 10 <= 10\leq 1010 and unfortunately underestimates that for Bo > 10. It remains an open problem to explain why, for large values of Bo, this inequality is not fulfilled.

CONCLUSIONS

The numerical experiments reported in this paper, as well as some others which are in progress and refer to FB problems from laminar boundary layer theory, entitle us to think that the main difficulty with this type of problem consists in the high nonnormality of their Jacobian matrices. The collocation method implemented by b v p 4 c b v p 4 c bvp4cb v p 4 cbvp4c, due to its global character, overcomes to some extent these difficulties. As a general remark, we have to observe that the problems corresponding to microgravity conditions, i.e., Bond number closed to zero, are easier to be solved and the results are much more accurate than those obtained when this number is increased.

REFERENCES

  1. R. Finn, "Capillary Free Surfaces," Preprint 7/1981 (Institut für Angewante Mathematik Univ. Bonn, 1981)
  2. C.-I. Gheorghiu, Studia Univ. "Babes-Bolyai," Mathematica L , 61 L , 61 L,61\mathbf{L}, 61L,61 (2005).
  3. J. Kierzenka and L. F. Shampine, ACM Trans. Math. Softw. 27, 299 (2001).
  4. A. Miele, A. K. Aggarwal and J. L. Tietze, J. Comput. Phys. 15, 117 (1974).
  5. L. N. Trefethen and D. Bau III, Numerical Linear Algebra (SIAM, 1997).

  1. 1 1 ^(1){ }^{1}1 The text was submitted by the author in English.
2008

Related Posts