3. :, c, ν. 4. Burgers : t + c x = ν 2 u x 2, (3), ν. 5. : t + u x = ν 2 u x 2, (4), c. 2 u t 2 = c2 2 u x 2, (5) (1) (4), (1 Navier Stokes,., ν. t +

B: 2016 12 2, 9, 16, 2017 1 6 1,.,,,,.,.,,,., 1,. 1. :, ν. 2. : t = ν 2 u x 2, (1), c. t + c x = 0, (2). e-mail: iwayama@kobe-u.ac.jp,. 1

3. :, c, ν. 4. Burgers : t + c x = ν 2 u x 2, (3), ν. 5. : t + u x = ν 2 u x 2, (4), c. 2 u t 2 = c2 2 u x 2, (5) (1) (4), (1 Navier Stokes,., ν. t + u x = 1 p ρ x + ν 2 u x 2, (6) 2 gnuplot png,,,., gnuplot png( ),.,. 2.1 png gnuplot, gnuplot. gnuplot 2

, gnuplot. 1.,, test.plt.,. set xrange[0:2*pi] set yrange[-1:1] set xlabel x set ylabel u(x,t) plot sin(x) title sin(x) set term png set output sincurve.png replot 2. gnuplot,, sin(x), sincurve.png, ( 1 ). gnuplot> cd test.plt gnuplot> load test.plt 2.2 gif 1. test2.plt), gunplot. * 1 set xrange[0:2*pi] set yrange[-1:1] set xlabel x set ylabel u(x,t) set term png set size 1, 1 set output sincurve_000.png *1 eps,, png. 3

1 sin(x) 0.5 u(x,t) 0-0.5-1 0 1 2 3 4 5 6 x 1 sin(x). plot sin(x) title t=0 set output sincurve_001.png plot sin(x+pi/8.) title t=1 set output sincurve_002.png plot sin(x+2.*pi/8.) title t=2... set output sincurve_015.png plot sin(x+15.*pi/8.) title t=15 2. ImageMagic convert, png gif. * 2 > convert sincurve_???.png sincurve.gif *2? 1. 4

3. sincurve.gif, gif. 2.3 gif FORTRAN, gnuplot, png., convert, png. 1. FORTRAN (sample.f). * 3!! sample.f! ( )! : 2015.06.23, 2016.12.2! 0 \le x \le Lx N! 0 \le t \le t_max dt j_max=t_max/dt! t_out (j_out=t_out/dt )! c :! pi:! u(x,t)=sin((x-x_0)+c(t-t_0))!--------------------------------------------------------------------- implicit none integer N, i, j, k, j_max, j_out real*8 t_max, dt, t_out, t_0 real*8 dx, x_0 real*8 pi, Lx, c parameter (pi=3.141592, Lx=2.D+0*pi) parameter (N=256, t_max=40.0d0, t_out=.25d+0, dt=1.d-2) *3, 0 x L x N, 0 t t max dt sin((x x 0 ) + c(t t 0 )) t out (fort.???). L x = 2π, N = 256, c = 1, t max = 40, dt = 10 2, t out = 0.25, x 0 = 0, t 0 = 0., fort.100 fort.260 161. 5

real*8 real*8 x(0:n+1), t(0:int(t_max/dt)+1) u(0:n+1) dx=lx/dble(n) c=1.d+0 x_0=0.d+0 t_0=0.d+0 j_max=int(t_max/dt) j_out=int(t_out/dt)! k=0 j=0 t(j)=dble(j)*dt do i=1,n x(i)=dble(i)*dx u(i)=dsin((x(i)-x_0)+c*(t(j)-t_0)) enddo! do i=1,n write(k+100, (3f10.3) ) x(i), t(j), u(i) enddo! do j=1, j_max do i=1, n x(i)=dble(i)*dx t(j)=dble(j)*dt u(i)=dsin((x(i)-x_0)+c*(t(j)-t_0)) enddo! 6

if (mod(j, j_out)==0) then k=k+1 do i=1,n write(k+100, (3f10.3) ) x(i), t(j), u(i) enddo endif end do stop end 2. sample.f, fort.100 fort.260 161. > gfortran sample.f >./a.out 3. gnuplot (sample.plt) *4, fort.100 fort.260. (mov100.png mov260.png.) gnuplot> load sample.plt sample.plt. 6, gnuplot FORTRAN DO. set xrange [0:2*pi] set yrange [-1:1] set term png set size square set xlabel x *4,. /home3/iwayama/exp_16/pde/exp1/sample.plt 7

set ylabel u(x) do for [i = 100:260] { <--- ii = sprintf("%03d", i) iii=sprintf("%03d", i-100) outfile = "mov". ii. ".png" infile = "fort.". ii set output outfile plot infile u 1:3 w l t "t=".iii } 4. convert, mov100.png mov260.png move.gif. > convert mov???.png mov.gif 2.4 sample.f diffusion.f,, 1, { 1 4πν(t + t0 ) exp (x x 0) 2 } 4ν(t + t 0 ). ν x 0, t 0, dt, t out,. sample.plt diffusion.plt, png gif. gif e-mail iwayama@kobe-u.ac.jp. (7) 3 1, t = ν 2 u x 2, (8) 8

., ν. 3.1 : (8), 0 x L x. N., x i, x x i = i x, x = L x N, (9).,, j, t t j = j t, (10). (8), x, t,, (8), i, j : : u(x, t) u(x i, t j ) (11) (8),.. x i 2 u/ x 2 x i 1, x i, x i+1 3 u. x, u(x i±1, t j ) Taylor u(x i±1, t j ) = u(x i, t j ) ± x x + 1 2 2 u x 2 ( x)2 + O( x 3 ), (12)., O( x 3 ) ( x) 3., 2 u/ x 2 2 u x 2 = u(x i+1, t j ) + u(x i 1, t j ) 2u(x i, t j ) ( x) 2 + O( x 2 ) (13). *5, ( x) 2. 2 u x 2 = u(x i+1, t j ) + u(x i 1, t j ) 2u(x i, t j ) ( x) 2 (14) *5 O( x), O( x) 2. 9

:, Euler., t 1, t = u(x i, t j+1 ) u(x i, t j ) t (15). :, (8), 2, 1 u(x i, t j+1 ) u(x i, t j ) t = ν u(x i+1, t j ) + u(x i 1, t j ) 2u(x i, t j ) ( x) 2, (16), u(x i, t j+1 ) = u(x i, t j ) + ν t ( x) 2 {u(x i+1, t j ) + u(x i 1, t j ) 2u(x i, t j )},(17), t j u, t j+1 u. 3.2 (17),.,., 0 x L x., u(0, t) = u(l x, t), (18). * 6 Gauss, u(x, 0) = exp [ a ( x L ) ] 2 x 2 (19). *6 u(x 0, t j ) = u(x N, t j ). 10

ν = 2 10 3, L x = 2π, N = 256, t = 1.0 10 2, a = 10. (diffusion sample.f).. gnuplot, diffusion sample.plt. * 7,., diffusion theory.f. * 8 3.3 : diffusion sample.f c1 c Euler c : Takahiro IWAYAMA c 2014.07.27 c 2015.07.07 c 2016.12.06--08 c, x, t 2 c. c u(x_i, t_j) -> u(i). u(i) u. c---------------------------------------------- c23456789012345678901234567890123456789012345678901234567890123456789012 program diffusion_sample implicit none integer i, j, k, N real*8 t, dt, t_max, t_out *7 /home3/iwayama/exp_16/pde/exp2/. *8 fort.500. gunplot u. fort.150 fort.550., (19), < x <.,. 11

real*8 x, dx real*8 pi, nu, Lx, a parameter (nu=2.0d-3, a=10.d+0) parameter (N=256) real*8 u(0:n+1)!. real*8 d2u(1:n)!u x 2. c parameters for output integer t_step, j_out c t_max=100.0d0!t_max dt=1.0d-2! t_k=k*dt t_out=1.0d+0!t_out t_step=int(t_max/dt) j_out=int(t_out/dt) pi=acos(-1.0d+0) Lx=2.0D+0*pi dx=lx/dble(n)!!! c k=0 t=dble(0) do i=1, N x=dble(i)*dx u(i)=dexp(-a*(x-lx/dble(2))**2) enddo c call bound_cond(u,n) c do i=0, N 12

enddo write(k+100,100) real(i)*dx, t, u(i) c ---------------------------------- do j=1, t_step t=t+dt c u 2 call second_deriv(u, d2u, dx, N) c do i=1, N c Euler!<-- enddo c call bound_cond(u, N) c if (mod(j, j_out)==0) then k=k+1 do i=0, N write(k+100,100) real(i)*dx, t, u(i) enddo close(k+100) endif enddo 100 format (3(1x,e12.5)) stop end c---------------------------------------------------------------------- subroutine second_deriv(u,d2u,dx,n) c u x 2 integer N, i 13

real*8 u(0:n+1)! u real*8 d2u(1:n)! u 2 real*8 dx do i=1,n end do!<-- return end c----------------------------------------------------------------------- subroutine bound_cond(u,n) c integer N real*8 u(0:n+1) u(0)=u(n) u(n+1)=u(1) return end c------------------------------------------------------------------------ 3.4 ν,.,. gif.,, L x, N, t).. (ν. von Neumann, 14

.) 3.5 von Neumann t, t,., ν x, t, : t ( x)2 2ν. (20) (20) von Neumann. (20), t, ν x. * 9 *9 u(x i, t j+1 ) = λu(x i, t j )e ikx j x Fourier, λ, Euler λ Euler = 1 + 2 ν t (cos(k x) 1) ( x) 2 = 1 4 ν t k x sin2 ( x) 2 2 (21). von Neumann 1 λ 1 (22), (21) 2 ν > 0., ν t ( x) 2 1 2 (23)., u(x, t + t) = λu(x, t)e ikx, λ theor = e k2 ν t (24). (21) (24), (21), k x 1, λ Euler = 1 νk 2 t + O((k x) 2 ) (25). (24), k 2 ν t 1, λ theor = 1 νk 2 t + O((k 2 ν t) 2 ) (26).. 15

4 *10 1,., c. t = c x, (27) 4.1 (27) x ct f(x ct) : u(x, t) = f(x ct). (28) (28) c x c, c x c. 4.2 : (27), 0 x L x. N., x i, x x i = i x, x = L x N, (29).,, j, t t j = j t, (30) *10, /home3/iwayama/exp_16/pde/exp3/. > cp /home3/iwayama/exp_16/pde/exp3/a*.. 16

. (27), x, t,, (27), i, j : u(x, t) u(x i, t j ) (31) : (8),.. x i / x x i 1, x i, x i+1 3 u. x, u(x i±1, t j ) Taylor u(x i±1, t j ) = u(x i, t j ) ± x x + 1 2 2 u x 2 ( x)2 + O( x 3 ), (32)., O( x) 3 ( x) 3., / x, 3 : 1. 2. 3. x = u(x i+1, t j ) u(x i, t j ) + O( x). (33) x x = u(x i, t j ) u(x i 1, t j ) + O( x). (34) x x = u(x i+1, t j ) u(x i 1, t j ) + O( x 2 ). (35) 2 x,, 3.,, O( x) ( x 1 ), O( x 2 ) ( x 2 ). :, Euler., t 1,. t = u(x i, t j+1 ) u(x i, t j ) t 17 (36)

:, (27), 1, 1 1. 2. u(x i, t j+1 ) = u(x i, t j ) c t x {u(x i+1, t j ) u(x i, t j )}. (37) u(x i, t j+1 ) = u(x i, t j ) c t x {u(x i, t j ) u(x i 1, t j )}. (38) 2, 1 3. u(x i, t j+1 ) = u(x i, t j ) c t 2 x {u(x i+1, t j ) u(x i 1, t j )}. (39), t j u, t j+1 u. 4.3 (37) (39),.,.,. * 11. L x /2, 1, u(0, t) = u(l x, t), (40) ( ) 4πx u(x, 0) = sin, (41) c = 0.1, L x = 2π, N = 256, t = 1.0 10 2. L x *11 u(x 0, t j ) = u(x N, t j ). 18

(advection.f), (advection sample.f).. gnuplot, advection.plt * 12.,., advection theory.f. gif.,,,, t). :,.. ( CFL, von Neumnann (42).). ( ), t = c x + ν 2 u x 2 (42),, ν = c x/2,. (42) 2. Courant-Friedrichs-Lewy(CFL). * 13 *12,,,. *13 CFL,. (, CFL.), von Neumann. u(x l, t j ) = λ j û 0 e ikx l, x l = l x, t j = j t (37), (38), (39). ( ) λ forward = 1 γ e ik x 1 = 1 + γ γe ik x, (43) ( λ backward = 1 γ 1 e ik x) = 1 γ + γe ik x, (44) λ central = 1 iγ sin k x, (45). γ c t/ x Courant. 2 λ forward 2 = 1 + 4γ(1 + γ) sin 2 k x, (46) λ backward 2 = 1 4γ(1 γ) sin 2 k x (47) λ central 2 = 1 + γ 2 sin 2 k x, (48) 19

5 *14 1, 2 u t 2., c * 15. = c2 2 u x 2, (50) 5.1 (50) x ct x + ct f(x ct), g(x + ct) : u(x, t) = f(x ct) + g(x + ct). (51) (51) d Alembert, c f(x ct) x c, g(x + ct) x c. f, g,. 5.2 (50) (8), (50) 2, (8) 1., 2 1. (46), (48) γ > 0 λ > 1., (47) 0 < γ 1 (49) von Neumann. (49) CFL.,, γ,. *14, /home3/iwayama/exp_16/pde/exp4/. > cp /home3/iwayama/exp_16/pde/exp3/m*.. *15. 20

., (50) 1,. v(x, t) t., (50) t = v, v t = c2 2 u x 2, (52a) (52b),,, Euler. 2, 1, u(x i, t j+1 ) u(x i, t j ) = v(x i, t j ), t (53a) v(x i, t j+1 ) v(x i, t j ) = c 2 u(x i+1, t j ) + u(x i 1, t j ) 2u(x i, t j ) t ( x) 2, (53b) u(x i, t j+1 ) = u(x i, t j ) + v(x i, t j ) t v(x i, t j+1 ) = v(x i, t j ) + c 2 t ( x) 2 {u(x i+1, t j ) + u(x i 1, t j ) 2u(x i, t j )}, (54a) (54b), t j u, v, t j+1 u, v., u, v. * 16 5.3 (54),.,. 0 x L x., (x, t) x = 0, (x = 0, L x ) (55). * 17 *16 2, 2. *17 u(x 0, t j ) = u(x 1, t j ), u(x N+1, t j ) = u(x N, t j ). 21

Gauss, u(x, 0) = exp [ a ( x L ) ] 2 x, v(x, 0) = 0, (56) 2. c = 1 10 1, L x = 2π, N = 256, t = 1.0 10 2, a = 10, 0 t 150. (wave sample.f).... gnuplot, wave.plt. * 18 5.4 1.,, u(0, t) = 0, u(l x, t) = 0 (57) (50),. 2. Euler,. Adams- Bashforth, ( ). 3. ( )Euler, von Neumann. * 19 5.5 1.,. *18 /home3/iwayama/exp_16/pde/exp4/. *19 : von Neumann 3.5. 22

c1 c Euler c : Takahiro IWAYAMA c 2014.07.27 c 2015.07.07 c 2016.12.06--08 c---------------------------------------------- c23456789012345678901234567890123456789012345678901234567890123456789012 program diffusion implicit none integer i, j, k, N real*8 t, dt, t_max, t_out real*8 x, dx real*8 pi, nu, Lx, a! <- nu -> c parameter (nu=2.0d-3, a=10.d+0)! <- c parameter (N=256) real*8 u(0:n+1)!! <- v real*8 d2u(1:n)!u x 2 c parameters for output integer t_step, j_out c t_max=100.0d0!t_max! <- dt=1.0d-2! t_k=k*dt t_out=1.0d+0!t_out t_step=int(t_max/dt) j_out=int(t_out/dt) pi=acos(-1.0d+0)! 23

Lx=2.0D+0*pi dx=lx/dble(n)!! c k=0 t=dble(0) do i=1, N x=dble(i)*dx u(i)=dexp(-a*(x-lx/dble(2))**2)! <- v enddo c call bound_cond(u,n) c do i=0, N write(k+100,100) real(i)*dx, t, u(i) enddo c ---------------------------------- do j=1, t_step t=t+dt c u 2 call second_deriv(u, d2u, dx, N) c do i=1, N c Euler u(i)=u(i)+nu*d2u(i)*dt! <-u, v enddo c call bound_cond(u, N) c if (mod(j, j_out)==0) then k=k+1 do i=0, N 24

write(k+100,100) real(i)*dx, t, u(i) enddo close(k+100) endif enddo 100 format (3(1x,e12.5)) stop end c---------------------------------------------------------------------- subroutine second_deriv(u,d2u,dx,n) c u x 2 integer N, i real*8 u(0:n+1) real*8 d2u(1:n) real*8 dx do i=1,n d2u(i)=(u(i+1)+u(i-1)-dble(2)*u(i))/dx**2 end do return end c----------------------------------------------------------------------- subroutine bound_cond(u,n) c integer N real*8 u(0:n+1) u(0)=u(n) u(n+1)=u(1)! <-! <- 25

return end c------------------------------------------------------------------------,, 211,, 2006.,, 6,, 1989. 26