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

B:,, 2017 12 1, 8, 15, 22 1,.,,,,.,.,,,., 1,. 1. :, ν. 2. : u t = ν 2 u x 2, (1), c. u t + c u x = 0, (2), ( ). 1

3. :, c, ν. 4. Burgers : u t + c u x = ν 2 u x 2, (3), ν. 5. : u t + u u x = ν 2 u x 2, (4), c. 2 u t 2 = c2 2 u x 2, (5) (1) (4), (1 Navier Stokes,., ν. u t + u 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.f90). * 3!! sample.f90! ( )! : 2015.06.23, 2016.12.2, 2017.11.30! 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 :: i, j, k, j_max, j_out real(8) :: t_0 real(8) :: dx, x_0 real(8) :: c integer, parameter :: N=256 real(8), parameter ::pi=3.141592, Lx=2.D+0*pi *3, 0 x L x N, 0 t t max dt sin((x x 0 ) + c(t t 0 )) t out (fort????.dat). L x = 2π, N = 256, c = 1, t max = 40, dt = 10 2, t out = 0.25, x 0 = 0, t 0 = 0., fort 0000.dat fort 0160.dat 161. 5

real(8), parameter :: t_max=40.0d0, t_out=.25d+0, dt=1.d-2 real(8) :: x(0:n+1), t, 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=t_0 do i=0,n x(i)=dble(i)*dx u(i)=dsin((x(i)-x_0)+c*(t-t_0)) enddo! call save_data(k, x, t, u, n)! do j=1, j_max t=t+dt do i=0, n x(i)=dble(i)*dx u(i)=dsin((x(i)-x_0)+c*(t-t_0)) enddo! if (mod(j, j_out)==0) then 6

k=k+1 call save_data(k, x, t, u, n) endif end do stop end program!---------------------------------------------------------------- subroutine counter(i, data_number) integer :: i, o_4, o_3, o_2, o_1 character(len=10) :: i_data= 0123456789 character(4) :: data_number o_4=i/1000 o_3=(i-o_4*1000)/100 o_2=(i-o_4*1000-o_3*100)/10 o_1= i-o_4*1000-o_3*100-o_2*10 data_number=i_data(o_4+1:o_4+1)// & i_data(o_3+1:o_3+1)// & i_data(o_2+1:o_2+1)// & i_data(o_1+1:o_1+1) return end subroutine counter!---------------------------------------------------------------------------- subroutine save_data(k, x, t, u, n) integer :: k, i, n real(8) :: x(0:n), t, u(0:n) character(4) :: data_number 7

call counter(k, data_number) do i=0, n open(50,file= fort_ //data_number//.dat, & status= unknown ) write(50, *) x(i), t, u(i) enddo end subroutine save_data 2. sample.f90, fort 0000.dat fort 0160.dat 161. > gfortran sample.f90 >./a.out 3. gnuplot (sample.plt) *4, fort 0000.dat fort 0160.dat. (mov 0000.png mov 0160.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_17/pde/exp1/sample.plt 8

set ylabel u(x) do for [i = 0:160] { <--- ii = sprintf("%04d", i) outfile = "mov_". ii. ".png" infile = "fort_". ii. ".dat" set output outfile plot infile u 1:3 w l t "t=".ii } 4. convert, mov 0000.png mov 0160.png move.gif. > convert mov_????.png mov.gif 2.4 sample.f90 diffusion.f90,, 1, { 1 4πν(t + t0 ) exp (x x 0) 2 } 4ν(t + t 0 ) (7). ν x 0, t 0, dt, t out,. sample.plt diffusion.plt, png gif., t 0,. * 5 gif e-mail iwayama@kobe-u.ac.jp. *5 t 0 = 0 t = 0 (7).. 9

3 1,., ν. u t = ν 2 u x 2, (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 ) ± u 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) 10

. *6, ( 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), Euler., t 1, u 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). * 7 *6 O( x), O( x) 2. *7 u(x 0, t j ) = u(x N, t j ). 11

Gauss, u(x, 0) = exp [ a ( x L ) ] 2 x 2 (19). ν = 2 10 3, L x = 2π, N = 256, t = 1.0 10 2, a = 10. (diffusion sample.f90).. gnuplot, diffusion sample.plt. * 8,., diffusion theory.f90. * 9 3.3 : diffusion sample.f90!1! Euler! : Takahiro IWAYAMA! 2014.07.27! 2015.07.07! 2016.12.06--08! 2017.12.2!---------------------------------------------- program diffusion implicit none *8 /home3/iwayama/exp_17/pde/exp2/. *9 theory 0000.dat. gunplot u., fort 0150.dat theory 0150.dat., (19), < x <.,. 12

integer :: i, j, k real(8) :: t, dt, t_max, t_out real(8) :: dx real(8) :: pi, Lx, nu integer, parameter :: N=256 real(8), parameter :: c=1.0d-1, a=10.d+0 real(8) :: x(0:n+1)!x real(8) :: u(0:n+1)! real(8) :: v(0:n+1)! real(8) :: d2u(1:n)!u x 2! parameters for output integer :: t_step, j_out! 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) nu=2.0d-3! pi=acos(-1.0d+0) Lx=2.0D+0*pi dx=lx/dble(n)!!!! k=0 t=dble(0) do i=1, N x(i)=dble(i)*dx u(i)=dexp(-a*(x(i)-lx/dble(2))**2) enddo 13

! call bound_cond(u,n)! call save_data(k, x, t, u, n)! ---------------------------------- do j=1, t_step t=t+dt! u 2 call second_deriv(u, d2u, dx, N)! do i=1, N! Euler enddo! call bound_cond(u, N)! if (mod(j, j_out)==0) then k=k+1 call save_data(k, x, t, u, n) endif enddo!<-- stop end program diffusion!---------------------------------------------------------------------- subroutine second_deriv(u, d2u, dx, N)! u x 2 integer :: N, i real(8) :: u(0:n+1) 14

real(8) :: real(8) :: d2u(1:n) dx do i=1,n d2u(i)=(u(i+1)+u(i-1)-dble(2)*u(i))/dx**2 end do return end subroutine second_deriv!----------------------------------------------------------------------- subroutine bound_cond(u,n)! integer :: N real(8) :: u(0:n+1) u(0)=u(n) u(n+1)=u(1) return end subroutine bound_cond!---------------------------------------------------------------- subroutine counter(i, data_number) integer :: i, o_4, o_3, o_2, o_1 character(len=10) :: i_data= 0123456789 character(4) :: data_number o_4=i/1000 o_3=(i-o_4*1000)/100 o_2=(i-o_4*1000-o_3*100)/10 o_1= i-o_4*1000-o_3*100-o_2*10 data_number=i_data(o_4+1:o_4+1)// & i_data(o_3+1:o_3+1)// & i_data(o_2+1:o_2+1)// & 15

i_data(o_1+1:o_1+1) return end subroutine counter!---------------------------------------------------------------------------- subroutine save_data(k, x, t, u, n) integer :: k, i, n real(8) :: x(0:n+1), t, u(0:n+1) character(4) :: data_number call counter(k, data_number) do i=1, n open(50,file= fort_ //data_number//.dat, & status= unknown ) write(50, *) x(i), t, u(i) enddo end subroutine save_data 3.4 ν,.,. gif.,, L x, N, t).. (ν. von Neumann,.) 16

3.5 von Neumann t, t,., ν x, t, : t ( x)2 2ν. (20) (20) von Neumann. (20), t, ν x. * 10 *10 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).. 17

4 *11 1,., c. u t = c u 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) *11, /home3/iwayama/exp_16/pde/exp3/. > cp /home3/iwayama/exp_16/pde/exp3/a*.. 18

. (27), x, t,, (27), i, j : u(x, t) u(x i, t j ) (31) : (8),.. x i u/ 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 ) ± u x x + 1 2 2 u x 2 ( x)2 + O( x 3 ), (32)., O( x) 3 ( x) 3., u/ x, 3 : 1. 2. 3. u x = u(x i+1, t j ) u(x i, t j ) + O( x). (33) x u x = u(x i, t j ) u(x i 1, t j ) + O( x). (34) x u 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,. u t = u(x i, t j+1 ) u(x i, t j ) t 19 (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),.,.,. * 12. 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 *12 u(x 0, t j ) = u(x N, t j ). 20

(advection.f), (advection sample.f).. gnuplot, advection.plt * 13.,., advection theory.f. gif.,,,, t). :,.. ( CFL, von Neumnann (42).). ( ), u t = c u x + ν 2 u x 2 (42),, ν = c x/2,. (42) 2. Courant-Friedrichs-Lewy(CFL). * 14 *13,,,. *14 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) 21

5 *15 1, 2 u t 2., c * 16. = 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.,, γ,. *15, /home3/iwayama/exp_16/pde/exp4/. > cp /home3/iwayama/exp_16/pde/exp3/m*.. *16. 22

., (50) 1,. v(x, t) u t., (50) u 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. * 17 5.3 (54),.,. 0 x L x., u(x, t) x = 0, (x = 0, L x ) (55). * 18 *17 2, 2. *18 u(x 0, t j ) = u(x 1, t j ), u(x N+1, t j ) = u(x N, t j ). 23

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. * 19 5.4 1.,, u(0, t) = 0, u(l x, t) = 0 (57) (50),. 2. Euler,. Adams- Bashforth, ( ). 3. ( )Euler, von Neumann. * 20 5.5 1.,. *19 /home3/iwayama/exp_16/pde/exp4/. *20 : von Neumann 3.5. 24

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)! 25

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 26

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)! <-! <- 27

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