11 Poisson kanenko@mbkniftycom alexeikanenko@docomonejp http://wwwkanenkocom/
, u = f, ( u = u+f u t, u = f t ) 1 D R 2 L 2 (D) := {f(x,y) f(x,y) 2 dxdy < )} D D f,g L 2 (D) (f,g) := f(x,y)g(x,y)dxdy (L 2 ) D f := (f,f) f L 2 ( ) L 2 (D) ( ) L 2 (D), l 2 = {(x 1,x 2,);x 2 1 +x2 2 + < }
Poisson - 1 (numb-1f) 2 D Dirichlet, ϕ j = λ j ϕ j, (D ), ϕ j D = 0 λ j ϕ j,,, {ϕ j } j=1 ie (ϕ i,ϕ j ) = δ ij, f(x,y) = c j ϕ j (x,y), c j = (f,ϕ j ) := f(x,y)ϕ j (x,y)dxdy j=1 D Fourier Fourier [0,1] [0,1] Poisson ( [a 0,a 1 ] [b 0,b 1 ] OK) (sinmπxsinnπy) = (m 2 +n 2 )π 2 sinmπxsinnπy, {sinmπxsinnπy} m,n=1, (m2 +n 2 )π 2 Q o 2sinmπxcosnπx 1 1 sin 2 1 cos2mπx mπxdx = dx = 1 0 0 2 2 1 { 2sinmπx} m=1
u(x,y) = m,n=1 u mn sinmπxsinnπy, f(x,y)= f mn sinmπxsinnπy, f mn = 4 m,n=1 D 3 f(x,y)sinmπxsinnπydxdy (m 2 +n 2 )π 2 u mn sinmπxsinnπy = f mn sinmπxsinnπy m,n=1 u mn = u(x,y) = m,n=1 f mn (m 2 +n 2 )π 2 m,n=1 f mn (m 2 +n 2 )π 2 sinmπxsinnπy,, 1 1, f(x,y)sinmπxsinnπydxdy 0 0
Poisson - 2 (numb-2f) [a 0,a 1 ] [b 0,b 1 ] Poisson ( 2 u x 2 + 2 u ) y 2 = f u(x+h x,y)+u(x h x,y) 2u(x,y) h 2 x u(x,y +h y)+u(x,y h y ) 2u(x,y) = f(x,y) h 2 y h x = a 1 a 0 M, h y = b 1 b 0 N, f ij = f(a 0 +ih x,b 0 +jh y ), u ij = u(a 0 +ih x,b 0 +jh y ), 2 + 2 1 0 1 0 0 h 2 x h 2 y h 2 x h 2 u y 11 f 11 1 2 + 2 1 h 2 x h 2 x h 2 y h 2 x u 21 f 21 0 1 h u M 1,1 f M 1,1 2 y 0 0 u 12 = f 12 1 h 2 y 0 u M 1,2 f M 1,2 1 h 2 x 0 0 1 0 1 2 + 2 h 2 y h 2 x h 2 x h 2 u y M 1,N 1 f M 1,N 1 4
Poisson - 3 (numb-3f) Poisson u = f D L 2 (D) {ϕ n } n=1 ( u, ϕ j ) = ( u,ϕ j ) = (f,ϕ j ) (, (f,g) = f(x,y)g(x,y)dxdy L 2 D ) u = u i ϕ i i=1 ( ϕ i, ϕ j )u i = f j, j = 1,2, i=1 u i Galerkin ( ϕ i ) {ϕ i } N i=1 u N ϕ i (Finite Element Method, FEM) A = (( ϕ i, ϕ j )) N i,j=1 5
D 6 1 (ie x,y 1 ) P1 1 L 2 D D = M i=1 T i P 1,,P N P i T i1,,t imi P i1,,p imi 1 ϕ i i 1 P i1,,p imi 0 0 i 1 P i m i T i m P i i T Pi 1 i 1 T i 2 Pi 2 P i m i P i 1 P i P i 2 ϕ i ( ) 0 A ( ϕ i, ϕ j )
7 (numb-3f ) 9 10 11 12 5 6 7 8 ϕ 6 (x,y) 1 2 3 4 Delaunay (FreeFEM ) FreeFEM ~kanenko/numcal/freefem examples, hogeedp ~kanenko/numcal/freefem/freefem++ hogeedp
(Conjugate gradient method, CG ) Ax = b f(x) := (Ax b,ax b) f (x), f(x+h) = (Ax+h b,ax+h b) = (Ax b,ax b)+(h,ax b)+(ax b,h)+o( h 2 ) = f(x)+2(ax b,h)+o( h 2 ), h ie Ax b = 0, f ( ), f f, (Ax,x) ( ), 8
9 0 x 0 1 r 0 = b Ax 0 p 0 = r 0 ( ) k = 0,1,2,,N 1 k -1 α k = (r k,r k ) (Ap k,p k ), x k+1 = x k +α k p k, r k+1 = r k α k Ap k (= b Ax k+1 ) r k+1 ε b k -2 β k = (r k+1,r k+1 ), p (r k,r k ) k+1 = r k+1 +β k p k 0 r 0,,r j = p 0,,p j, j = 0,,k r k+1 = 0 x k+1 r j r N = 0 A ( ), r k+1
10 j < k (r k,r j ) = 0, (Ap k,p j ) = 0 (p k,r k ) = (r k,r k ) = (p k,r k 1 ) (r k,p j ) = 0 (j < k), (p k,r j r j 1 ) = 0 (j k), k k = 0 k, k+1, α k (r k+1,r k ) = (r k α k Ap k,r k ) = (r k,r k ) α k (Ap k,r k ) = (r k,r k ) α k (Ap k,p k )+α k β k 1 (Ap k,p k 1 ) = 0 (r k+1,p k ) = (r k,p k ) α k (Ap k,p k ) = (r k,r k ) α k (Ap k,p k ) = 0 (p k+1,r k+1 ) = (r k+1,r k+1 )+β k (p k,r k+1 ) = (r k+1,r k+1 ) (p k+1,r k ) = (r k+1,r k )+β k (p k,r k ) = β k (r k,r k ) = (r k+1,r k+1 ) 1 (Ap k+1,p k ) = (p k+1,ap k ) = (p k+1, (r k r k+1 )) α k = 1 {(p α k+1,r k ) (p k+1,r k+1 )} = 0 k, j < k r j p 0,,p j, (r k+1,r j ) = (r k α k Ap k,r j ) = (r k,r j ) α k (Ap k,r j ) = 0 (Ap k+1,p j ) = (p k+1,ap j ) = (r k+1 +β k p k, 1 α j (r j r j+1 )) = β k α j (p k,r j+1 r j ) = 0 QED
(Choleski) A L A = LL T LU A = (a ij ), L = (s ij ) A = LL T i i a ii = s 2 ik, a ij = s ik s jk (i < j) k=1 k=1 a 11 = s 2 11, a 12 = s 11 s 21,, a 1n = s 11 s n1, a 22 = s 2 21 +s2 22, a 23 = s 21 s 31 +s 22 s 32,, a 2n = s 21 s n1 +s 22 s n2 a ij = a ji L 1 s 11 = a 11, s i1 = a i1, i = 2,,n s 11 2 s 22 = a 22 s 2 21, s i2 = a i2 s 21 s i1, i = 3,,n s 22 k s kk = a kk k 1 j=1 s2 kj, s ik = a ik k 1 j=1 s ijs kj s kk 11, i = k+1,,n a 11 a 12 a 1n s 11 0 0 s 11 s 21 s n1 a 21 a 22 a 2n s 21 s 22 0 s 22 s n2 0 a n1 a n2 a nn s n1 s n2 s nn 0 0 s nn Choleski A
12 u(x,y,t) = v(x,y)w(t) Dirichlet, 1 2 u c 2 t 2 = 2 u x 2 + 2 u y 2 1 c 2v(x,y)w (t) = v(x,y)w(t), v(x, y) = w (t) v(x, y) c 2 = λ ( w(t) ) v(x,y) = λv(x,y) on D, v(x,y) D = 0 w +c 2 λw = 0 w(t) = c 1 cosc λt+c 2 sinc λt c λ/2π D = [0,1] [0,1] sinmπxsinnπy (m 2 +n 2 )π 2 ( eig2dfemf)
13 (Householder) (Lanczos),,, Lanczos 1 u 1 2 α 1 = (Au 1,u 1 ), v 2 = Au 1 α 1 u 1, β 1 = v 2, u 2 = v 2 /β 1 k u k, u k 1 v k+1 = Au k β k 1 u k 1 α k u k, β k = v k+1, u k+1 = v k+1 /β k P = (u 1,,u n ) B = P 1 AP AP = PB α 1 β 1 0 0 β 1 α 2 β 2 A(u 1,,u n ) = (u 1,,u n ) 0 β 2 0 βn 1 0 0 β n 1 α n Au 1 = α 1 u 1 +β 1 u 2, Au 2 = β 1 u 1 +α 2 u 2 +β 2 u 3,, Au k = β k 1 u k 1 +α k u k +β k u k+1, u k+1
14 k = 1,2,,n λ α 1 β 1 0 0 β 1 λ α 2 β 2 p k (λ) := 0 β 2 0 βk 1 0 0 β k 1 λ α k p k (λ) = (λ α k )p k 1 (λ) β 2 k 1 p k 2(λ) p 0 (λ) = 1, p 1 (λ) = 0 k = 1, λ, p k (λ) (p n (λ),,p 1 (λ),p 0 (λ)) Sturm ie a < b a b = [a,b]
(condition number), κ(a) = A A 1 L 2 - Ax A = sup x x 0 x x, A,, A = A,, A 1 A A 1 = A 1 = A, max λ i 1 i N κ(a) = min λ i = A A 1 i N 1, (Ax,x) = 0, A 1 A 15
16 Ax = b, A, b A, b, x A < 1 A 1, x κ(a) ( A x 1 κ(a) A / A A + b ) b B B < 1, I B, Neumann (I B) 1 = I +B +B 2 + (, I B I ) (I B) 1 = I +B +B 2 + 1+ B + B 2 + 1 = 1 B, B = A 1 A, (I +A 1 A) 1 1 1 A 1 A 1 1 A 1 A
, Ax = b (A+ A)(x+ x) = b+ b, 17 Ax+(A+ A) x = b A 1 A 1 Ax+(I +A 1 A) x = A 1 b x, x = (I +A 1 A) 1 A 1 ( Ax+ b) x (I +A 1 A) 1 A 1 ( A x + b ) A 1 1 A 1 A ( x + b ) A A 1 1 A 1 A ( A x + b A A ), b A x, κ(a) x A A 1 ( A x 1 A 1 A A + b ) A x κ(a) ( A 1 κ(a) A / A A + b ) QED b, A (, )
18, i a ii > j i a ij, Jacobi, 1 (preconditioning) Jacobi, Gauss-Seidel,, 1, (successive over-relaxation, SOR, ) x k+1 = x k +R k, 0 < ω < 1, x k+1 (1 ω)x k +ωx k+1 = x k +ωr k, ( ) 1 < ω < 2,,,, over-relaxation ω Jacobi Gauss-Seidel A, Jacobi M J (A), Gauss-Seidel M GS (A), ρ(m GS (A)) < ρ(m J (A)) < 1 ρ(m GS (A)) > ρ(m J (A)) > 1 ρ(m GS (A)) = ρ(m J (A)) = 1, Gauss-Seidel ( )
LAPACK (Linear Algebra Package) 19 LINPACK FORTRAN77, SCALAPACK, F95 lapack95, C clapack, C++ lapack++, JAVA ~kanenko/numcal/lapack/src,, BLAS (Basic Linear Algebra Subprograms) LAPACK, BLAS,,,,
20 1 Poisson Fourier numb-1f N, 2 Poisson numb-2f h, 3 Poisson numb-3f h, 4 FreeFEM
Q o FreeFEM JL Lions 21 /home/isstaff/kanenko/numcal/freefem++/bin /home/isstaff/kanenko/numcal/freefem++/examples* cd examples++ FreeFem++-x11 demoedp FreeFem numb-1f poisson1f, numb-2f poisson2f, numb-3f poisson3f, eig2dfemf ( ) wave2dfemf ( )
22 111 [0,1] [0,1] R 2 Poisson u = 1 Dirichlet Fourier 112 [0,1] [0,1] R 2 Poisson u = 1 Dirichlet 113 [0,1] [0,1] R 2 N Dirichlet Poisson u = 1