Relaxation scheme of Besse t t n = n t, u n = u(t n ) (n = 0, 1,,...)., t u(t) = F (u(t)) (1). (1), u n+1 u n t = F (u n ) u n+1 = u n + tf (u n )., t

Size: px

Start display at page:

Download "Relaxation scheme of Besse t t n = n t, u n = u(t n ) (n = 0, 1,,...)., t u(t) = F (u(t)) (1). (1), u n+1 u n t = F (u n ) u n+1 = u n + tf (u n )., t"

かずまさこしの
5 years ago
Views:

1 RIMS relaxation sheme of Besse splitting method Scilab Scilab Google Scilab Scilab Mathieu Colin Mathieu Colin 1

2 Relaxation scheme of Besse t t n = n t, u n = u(t n ) (n = 0, 1,,...)., t u(t) = F (u(t)) (1). (1), u n+1 u n t = F (u n ) u n+1 = u n + tf (u n )., t = t n+1/ (1) u n+1 u n = F (u n+1/ ), u n+1/ = un+1 + u n t (Crank-Nicolson scheme). () u(t + h) u(t) h u(t ± h) = u(t) ± hu (t) + h u (t) + O(h 3 ) = u (t) + O(h),. () u(t + h) u(t h) h = u (t) + O(h ) u n+1 = u n+1/ u n, (3) u n+1/ tf (u n+1/ ) = u n (4). F, u n u n+1/, (4)., F F (u) = Au, (4) (I ta)u n+1/ = u n,, u n u n+1/, (3) u n+1.

3 , A, g t u = Au + g(u)u (5). Crank-Nicolson scheme,. u n+1 u n t = Au n+1/ + g(u n )u n+1/, u n+1/ = un+1 + u n (6), u n, (I ta tg(u n ))u n+1/ = u n, u n+1/. (6), g(u n ) t = t n g(u(t)), t = t n+1/. g(u) t = t n+1/, Crank-Nicolson scheme., (Relaxation scheme of Besse). ϕ, (5) { ϕ = g(u) (7) t u = Au + ϕu (8). (7) t = t n, (8) t = t n+1/ ϕ n+1/ + ϕ n 1/ = g(u n ), u n+1 u n = Au n+1/ + ϕ n+1/ u n+1/ t., ϕ 1/ = g(u 0 )., (ϕ n 1/, u n ) ϕ n+1/ = g(u n ) ϕ n 1/, (I ta tϕ n+1/ )u n+1/ = u n, ϕ n+1/ u n+1/. u n+1 = u n+1/ u n u n+1, (ϕ n+1/, u n+1 ). 3

4 Schrödinger i t u = xu + V (x)u g(u)u, 0 x < L, 0 < t T, u(t, 0) = u(t, L), 0 < t T, u(0, x) = u 0 (x), 0 x < L relaxation scheme of Besse., V (x), g(u) = c u p (c R, p > 0). [0, L) N. x = L/N. x j = j x (j = 0, 1,..., N 1). t, t n = n t (n = 0, 1,,...). u(t, x) u n = [u n 0, u n 1,..., u n ] T, u n j = u(t n, x j )., T. n, u n N. u = [u 0, u 1,..., u ] T, g(u) = [g(u 0 ), g(u 1 ),..., g(u )] T. Laplacian u(x ± h) = u(x) ± hu (x) + h u (x) ± h3 3! u(3) (x) + O(h 4 ) u(x + h) + u(x h) u(x) = u (x) + O(h ). h, x N D xx = 1 ( x) , A = D xx diag [V (x 0 ), V (x 1 ),..., V (x )]. 4

5 w n+1/ + w n 1/ = g(u n ), w 1/ = g(u 0 ) i un+1 u n t + Au n+1/ + w n+1/ u u+1/ = 0, u n+1/ = un+1 + u n., u n+1 u n = (u n+1/ u n ), W n+1/ = diag (w n+1/ ) w n+1/ = g(u n ) w n 1/, ( iin + ta + tw n+1/) u n+1/ = iu n., I N N., u n u n+1/. u n+1 = u n+1/ u n, u n+1. Scilab. N u, v, u, v = x u j v j, u = u, u 1/ j=0. {u n }, L., n = 0, 1,,..., u n = u 0.,,, []. []. [1] C. Besse, Schéma de relaxation pour l équation de Schrödinger non linéaire et les systèmes de Davey et Stewartson, C. R. Acad. Sci. Paris Sér. I Math. 36 (1998) [] C. Besse, A relaxation scheme for the nonlinear Schrödinger equation, SIAM J. Numer. Anal. 4 (004)

6 Solitary waves for nonlinear Schrödinger equations Schrödinger. i t u + xu + g(u)u = 0, (x, t) R g(u) = c 3 u p (p > 0, c 3 > 0) ω > 0, k R. u ω,k (t, x) = exp(ikx ik t) exp(iω t)φ ω (x kt), ( ) 1/p ω c3 φ ω (x) =, d 3 = d 3 cosh(pωx) p + 1 u ω,k (0, x) = (ω/d 3 ) 1/p exp(ikx){cosh(pωx)} 1/p. g(u) = c u p + c 3 u p (p > 0, c, c 3 R) ω > 0, k R φ ω (x) = (. u ω,k (t, x) = exp(ikx ik t) exp(iω t)φ ω (x kt), ω d + d + d 3ω cosh(pωx) u ω,k (0, x) = exp(ikx) ( ) 1/p, d = c p +, d c3 3 = p + 1 ω d + d + d 3ω cosh(pωx) ) 1/p. 1

7 05/5/11 09:54:03 C:/Documents and Settings/Administrator/My Documents/Kyoto053/B clear; // Nonlinear Schrodinger equation // iu_t+u_{xx}+c3 u ^{p}u=0, 0<x<L, 0<t<T // with periodic boundary conditions // Constants c3=0; p=1; L=10; N=56; T=10; NT=0; dt=10^(-3); ome1=4; k1=0; ome=3; k=-1; dx=*l/n; x=[0:n-1]'*dx; // Space variable Dxx=-*speye(N,N)+[spzeros(1,N);speye(N-1,N)]... +[spzeros(n,1),speye(n,n-1)]; Dxx(1,N)=1; Dxx(N,1)=1; Dxx=Dxx/dx/dx; // Second Derivative in x function y=g(u) // Nonlinear term y=c3*abs(u)^(*p) d3=sqrt(c3/(p+1)); // Solitary waves function y=f1(x) y=(ome1/d3)^(1/p)*exp(%i*k1*x).*(cosh(p*ome1*x))^(-1/p) function y=f(x) y=(ome/d3)^(1/p)*exp(%i*k*x).*(cosh(p*ome*x))^(-1/p) u0=f1(x-l+4)+f(x-l-4); // Initial data IA=*%i*speye(N,N)+dt*Dxx; timer(); // Relaxation scheme of Besse t=0; u=u0; w=g(u); for n=1:nt, while t<=n*t/nt, w=*g(u)-w; v=(ia+dt*sparse(diag(w)))\(*%i*u); u=*v-u; clear v; t=t+dt; end clf(); // Animation plot(x,real(u),'g'); plot(x,imag(u),'y'); plot(x,abs(u),'b'); plot(x,-abs(u),'b'); Mass(n)=sum(abs(u)^)*dx; Energy(n)=-real(u'*Dxx*u)*dx/-c3*sum(abs(u)^(*p+))*dx/(*p+); end Timer=timer() plot(x,abs(u0),'r'); xset("window",1) plot(mass); xset("window",) plot(energy); ErrorMass=max(Mass)-min(Mass) ErrorEnergy=max(Energy)-min(Energy) Page : 1

8 Splitting method for nonlinear Schrödinger equations, splitting method, t u = A(u) + B(u), u(0) = u 0 u(t) = S(t)u 0, t v = A(v) X(t) t w = B(w) Y (t). Z L (t) = Y (t)x(t), Z S (t) = Y (t/)x(t)y (t/)., T > 0, t > 0, n t T n N S(n t)u 0 Z L ( t) n u 0 C t, S(n t)u 0 Z S ( t) n u 0 C( t)., C T u 0. Z L Lie formula, Z S Strang formula., Y (t) Z S ( t) n = Y ( t/)z L ( t) n Y ( t/)., Schrödinger { i t u = (1) xu + V (x)u g(u)u (t, x) R u(0, x) = u 0 (x) x R., V (x), g, g(z) = g( z ) (z C). g(u) = c u p (c R, p > 0)., (1) { t v = i () xv (t, x) R v(0, x) = v 0 (x) x R (3). () {., (3) t w = i[g(w) V (x)]w (t, x) R w(0, x) = w 0 (x) x R v(t, x) = F 1 [e itξ Fv 0 ](x) w(t, x) = exp ( it{g(w 0 (x)) V (x)} ) w 0 (x) (4)., g, x R t w(t, x) = R(ww) = R{i[g(w(t, x)) V (x)] w(t, x) } = 0., (t, x) R w(t, x) = w 0 (x)., g(z) = g( z ), (t, x) R g(w(t, x)) = g(w 0 (x)), (4). Scilab. 1

9 [1] C. Besse, B. Bidégaray, and S. Descombes, Order estimates in time of splitting methods for nonlinear Schrödinger equation, SIAM J. Numer. Anal. 40 (00) [] H. Holden, K. H. Karlsen, K.-A. Lie, and N. H. Risebro, Splitting methods for partial differential equations with rough solutions, Analysis and MATLAB programs, EMS Series of Lectures in Mathematics, European Mathematical Society (EMS), Zürich, 010. [3] H. Holden, K. H. Karlsen, N. H. Risebro, and T. Tao, Operator splitting for the KdV equation, preprint, arxiv: [4] C. Lubich, On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations, Math. Comp. 77 (008)

10 Fourier f(x) L., k Z, ξ k = kπ/l, f(x) Fourier ˆf(ξ k ) = 1 L L 0 f(x)e iξ kx dx (1) f(x) =.. k= ˆf(ξ k )e iξ kx () N, [0, L) N. x = L/N. x j = j x, f j = f(x j ) (j Z)., f j+n = f j (j Z), f(x) N [f 0, f 1,..., f ]. (1) ˆf(ξ k ) 1 f(x j )e iξ kx j x = L ˆf k j=0. ˆf k = 1 N j=0 f j ω jk, ω = exp ( ) πi N., ω N = 1, ˆfk+N = ˆf k (k Z), ˆf(ξk ) N [ ˆf 0, ˆf 1,..., ˆf ] [ ˆf N/+1,..., ˆf 1, ˆf 0, ˆf 1,..., ˆf N/ ]., () N/ k= N/+1 ˆf k ω jk = ˆf k ω jk., j = 0, 1,..., N 1 ˆf k ω jk = 1 N l=0 f l ω (j l)k = 1 N l=0 f l ω (j l)k, j, l = 0, 1,..., N 1, j l N 1, l j, ω j l 1 ω (j l)k = 1 ω(j l)n 1 ω j l = 0. 1

11 , l = j, ω (j l)k = N, () f j = ˆf k ω jk = N/ k= N/+1 ˆf k ω jk (j Z)., (), f (x). N/ k= N/+1 f (x) = k= iξ k ˆf(ξk )e iξ kx N/ iξ k ˆfk ω jk = iξ k ˆfk ω jk + k=n/+1 iξ k N ˆfk ω jk Scilab, N u = [u 0, u 1,..., u ], v=fft(u) N v = [v 0, v 1,..., v ], v k = j=0 u j ω jk., N v = [v 0, v 1,..., v ], u=ifft(v) N u = [u 0, u 1,..., u ], u j = 1 N. fft ifft. v k ω jk xi = [ξ 0, ξ 1,..., ξ N/, ξ N/+1,..., ξ 1 ] = π [0, 1,..., N/, N/ + 1,..., 1] L, p( x )u(x) ifft(p(%i*xi).*fft(u))., u = [u(x 0 ), u(x 1 ),..., u(x )], %i, p(z) v = [v 0, v 1,..., v ], p(v) p(v) = [p(v 0 ), p(v 1 ),..., p(v )], v w, v. w. v. w = [v 0 w 0, v 1 w 1,..., v w ]

12 05/5/11 09:57:13 C:/Documents and Settings/Administrator/My Documents/Kyoto053/S clear; // Nonlinear Schrodinger equation // iu_t+u_{xx}+c_3 u ^{p}u=0, 0<x<L, 0<t<T // with periodic boundary conditions // Constants c3=0; p=1; L=10; N=56; T=10; NT=0; dt=10^(-4); ome1=4; k1=0; ome=3; k=-1; dx=*l/n; x=[0:n-1]*dx; // Space variable xi=[0:n/,-n/+1:-1]*%pi/l; // Fourier variable function y=g(u) // Nonlinear term y=c3*abs(u)^(*p) function v=dx(u) // Derivative in x v=ifft(%i*xi.*fft(u)) d3=sqrt(c3/(p+1)); // Solitary waves function y=f1(x) y=(ome1/d3)^(1/p)*exp(%i*k1*x).*(cosh(p*ome1*x))^(-1/p) function y=f(x) y=(ome/d3)^(1/p)*exp(%i*k*x).*(cosh(p*ome*x))^(-1/p) u0=f1(x-l+4)+f(x-l-4); // Initial data timer(); // Splitting method with Strang formula t=0; u=u0; for n=1:nt, u=exp(%i*dt/*g(u)).*u; while t<=n*t/nt, u=ifft(exp(-%i*dt*xi^).*fft(u)); u=exp(%i*dt*g(u)).*u; t=t+dt; end u=exp(-%i*dt/*g(u)).*u; clf(); // Animation plot(x,real(u),'g'); plot(x,imag(u),'y'); plot(x,abs(u),'b'); plot(x,-abs(u),'b'); Mass(n)=sum(abs(u)^)*dx; Energy(n)=sum(abs(Dx(u))^)*dx/-c3*sum(abs(u)^(*p+))*dx/(*p+); end Timer=timer() plot(x,abs(u0),'r'); xset("window",1) plot(mass); xset("window",) plot(energy); ErrorMass=max(Mass)-min(Mass) ErrorEnergy=max(Energy)-min(Energy) Page : 1

163 KdV KP Lax pair L, B L L L 1/2 W 1 LW = ( / x W t 1, t 2, t 3, ψ t n ψ/ t n = B nψ (KdV B n = L n/2 KP B n = L n KdV KP Lax W Lax τ KP L ψ τ τ Cha

163 KdV KP Lax pair L, B L L L 1/2 W 1 LW = ( / x W t 1, t 2, t 3, ψ t n ψ/ t n = B nψ (KdV B n = L n/2 KP B n = L n KdV KP Lax W Lax τ KP L ψ τ τ Cha 63 KdV KP Lax pair L, B L L L / W LW / x W t, t, t 3, ψ t n / B nψ KdV B n L n/ KP B n L n KdV KP Lax W Lax τ KP L ψ τ τ Chapter 7 An Introduction to the Sato Theory Masayui OIKAWA, Faculty of Engneering,