|
数値解析講座
|
| 科学技術者の技術の向上を支援することを目的として、基礎的な数値解析講座として Newton法、Runge-kutta 法、Euler法、Bisection法についての解法を解説するとともに、それらのソースプログラムを紹介した。また、流出解析手法として、広く活用されてきている貯留関数法についても、解説とソースプログラムを紹介した。 |
|
●数値解析法の説明（PDF) Newton法 Runge-kutta 法 Euler法 Bisection法 ●ソースプログラム（fortran言語） 方程式X ： f(x)=a+X**P （P=3、a=3)
●計算結果（貯留関数法） 方程式X ソースプログラム
c program
open(10,file='out.csv')
c 方程式 f(x)=a + x**p
a=-3.0
p=3.0
eps=1.e-5
call newton(a,p,eps)
call runge(a,p,eps)
call euler(a,p,eps)
call bisection(a,p,eps)
pause 'end'
stop
end
c-----------------------------------------------------
subroutine newton(a,p,eps)
funk(x) = a + x**p
funk2(x)=p*x**(p-1)
write(10,*) 'newton'
write(10,*) 'f(x),x'
x = 3.0 !初期値
do i=1,100
ax=funk(x)/funk2(x)
x=x-ax
if(abs(funk(x)).le.eps) go to 1000
write(10,9999) funk(x),x
9999 format(f10.4,',',f10.4)
enddo
1000 continue
write(*,*) 'newton ans= ',x
write(10,'(a6,f10.3)') 'ans= ,',x
return
end
c-----------------------------------------------------
subroutine runge(a,p,eps)
funk(x) = a + x**p
write(10,*) 'runge'
write(10,*) 'f(x),x'
x = 0.1 !初期値
h = 0.1 !刻み幅
do i=1, 100
a1 = h * funk(x)
a2 = h * funk(x + a1/2.)
a3 = h * funk(x + a2/2.)
a4 = h * funk(x + a3)
x = x - (a1 + 2.*a2 + 2.*a3 + a4)/6.
if(abs(funk(x)).le.eps) go to 1000
write(10,9999) funk(x),x
9999 format(f10.4,',',f10.4)
end do
1000 continue
write(*,*) 'runge ans= ',x
write(10,'(a6,f10.3)') 'ans= ,',x
return
end
c-----------------------------------------------------
subroutine euler(a,p,eps)
funk(x)=a+x**p
write(10,*) 'euler'
write(10,*) 'f(x),x'
x=0.1 !初期値
h=0.1 !刻み幅
do i=1,50
x=x-funk(x)*h
if(abs(funk(x)).le.eps) go to 1000
write(10,9999) funk(x),x
9999 format(f10.4,',',f10.4)
enddo
1000 continue
write(*,*) 'euler ans= ',x
write(10,'(a6,f10.3)') 'ans= ,',x
return
end
c-----------------------------------------------------
subroutine bisection(a,p,eps)
funk(x)=a+x**p
write(10,*) 'bisection'
write(10,*) 'f(x),x'
dh=5.0 !上限値
dl=-1.0 !下限値
x=3.0 !初期値
3000 continue
write(10,9999) funk(x),x
9999 format(f10.4,',',f10.4)
if(abs(funk(x)).le.eps) go to 1000
if(funk(x).gt.0.0) go to 2000
dl=x
x=(dh+dl)/2.0
go to 3000
2000 dh=x
x=(dh+dl)/2.0
go to 3000
c
1000 continue
write(*,*) 'bisection ans= ',x
write(10,'(a6,f10.3)') 'ans= ,',x
return
end
c program
dimension r(10)
data r/3.0, 8.0,16.0,26.0,14.2,3.0,0.0,0.0,0.0,0.0/ ! 時間雨量
open(10,file='out.csv') !出力ファイル
write(10,*) '方程式（貯留関数法）'
+ ,'f(sq2)=r-(sq1＋sq2)/2-k*sq2**p+k*sq1**p を解く (Q=sq*A/3.6)'
eps=1.e-5
ak=8.0
p=0.6
a=10.8
write(10,*) '○定数'
write(10,*) 'k=,',ak
write(10,*) 'p=,',p
write(10,*) 'A=,',a,',km2'
write(10,*) 'TL=,',' 0.0',',hr'
write(10,*) 'f1=,',' 0.0'
write(10,*) 'Rsa=,',' 0.0',',mm'
c
write(10,*) '○計算結果'
c
write(*,*) 'Newton'
write(10,*) 'Newton,r(mm/hr),sq(mm/hr),Q(m3/s)'
sq1=0.001
do i=1,10
call newton(ak,p,a,r(i),sq1,eps,x)
sq1=x
enddo
c
write(*,*) 'Runge'
write(10,*) 'Runge,r(mm/hr),sq(mm/hr),Q(m3/s)'
sq1=0.001
do i=1,10
call runge(ak,p,a,r(i),sq1,eps,x)
sq1=x
enddo
c
write(*,*) 'Euler'
write(10,*) 'Euler,r(mm/hr),sq(mm/hr),Q(m3/s)'
sq1=0.001
do i=1,10
call euler(ak,p,a,r(i),sq1,eps,x)
sq1=x
enddo
c
write(*,*) 'Bisection'
write(10,*) 'Bisection,r(mm/hr),sq(mm/hr),Q(m3/s)'
sq1=0.001
do i=1,10
call bisection(ak,p,a,r(i),sq1,eps,x)
sq1=x
enddo
pause 'end'
stop
end
c----------------------------------------------------------
subroutine newton(ak,p,a,r,sq1,eps,x)
funk(sq2)=r-(sq1+sq2)/2-ak*sq2**p+ak*sq1**p
funk2(sq2)=-ak*p*sq2**(p-1)-0.5 ! = f’(sq2)
x = 0.0001 !初期値
do i=1,100
xn=funk(x)/funk2(x)
x=x-xn
if(abs(funk(x)).le.eps) go to 1000
enddo
1000 continue
write(10,9000) r,x,x*a/3.6
9000 format(3(',',f10.3))
write(*,*) r,x,x*a/3.6
end
c----------------------------------------------------------
subroutine runge(ak,p,a,r,sq1,eps,x)
funk(sq2) = r-(sq1+sq2)/2-ak*sq2**p+ak*sq1**p
x = 100.0 !初期値
h = -0.1 !刻み幅
do i=1, 100
a1 = h * funk(x)
a2 = h * funk(x + a1/2.)
a3 = h * funk(x + a2/2.)
a4 = h * funk(x + a3)
x = x - (a1 + 2.*a2 + 2.*a3 + a4)/6.
if(abs(funk(x)).le.eps) go to 1000
end do
1000 continue
write(10,9000) r,x,x*a/3.6
9000 format(3(',',f10.3))
write(*,*) r,x,x*a/3.6
return
end
c----------------------------------------------------------
subroutine euler(ak,p,a,r,sq1,eps,x)
funk(sq2) = r-(sq1+sq2)/2-ak*sq2**p+ak*sq1**p
x=100.0 !初期値
h=-0.1 !刻み幅
do i=1,50
x=x-funk(x)*h
if(abs(funk(x)).le.eps) go to 1000
enddo
1000 continue
write(10,9000) r,x,x*a/3.6
9000 format(3(',',f10.3))
write(*,*) r,x,x*a/3.6
return
end
c----------------------------------------------------------
subroutine bisection(ak,p,a,r,sq1,eps,x)
funk(sq2) = r-(sq1+sq2)/2-ak*sq2**p+ak*sq1**p
dh=1000.0 !上限値
dl=0.0 !下限値
x=dh !初期値
3000 continue
if(abs(funk(x)).le.eps) go to 1000
if(funk(x).gt.0.0) go to 2000
dh=x
x=(dh+dl)/2.0
go to 3000
2000 dl=x
x=(dh+dl)/2.0
go to 3000
c
1000 continue
write(10,9000) r,x,x*a/3.6
9000 format(3(',',f10.3))
write(*,*) r,x,x*a/3.6
return
end
|
|
数値解析講座
|