program MetodaN
C
parameter (maxstp=10000)
dimension xk(0:maxstp), yk(0:maxstp)
C
external funu, funv
C
C========================================================================
C
C ..Opis programu MetodaN..
C
C Program rozwiazuje przykladowe rownanie rozniczkowe
C
C x**3 y'' = y, [y(1)=1, y'(1)=0]
C
C na odcinku [eps, 1] metod/a Numerova-Cowella ...
C
C Funkcje funu i funv (external) obliczają prawą strone rownania
C y''= U(x) +V(x) y (tak jak wymaga algorytm Numerova).
C
C Procedura y1Num wylicza wartosc y1=y(x0+h)
C z dokladnoscia O(h**5); procedura y2Num oblicza ka/zde nast/epne
C y(k+1) na podstawie y(k) i y(k-1). Procedura Numerow organizuje
C calosc rachunkow wykonujac numstep krokow na przedziale
C [x0, xn]. Wyniki obliczen znajduja sie w tablicach xk(0:numstep)
C oraz yk(0:numstep) i sa drukowane do pliku 'numout'.
C
C Znak kroku ca/lkowania jest wybierany automatycznie w zaleznosci
C od tego czy xn > x0 czy tez x0 > xn. W ostatnim przypadku
C h < 0 (tak jak w przykladzie).
C
C Trzecia kolumna wynikow w pliku 'numout' zawiera iloraz
C y(x)/(x**(3/4)*exp(1./sqrt(x)))
C
C========================================================================
C
C ..actual number of steps..
numstp = 10000
if (numstp.gt.maxstp) stop 'NUM: Actual number of steps too big!'
C
C ..initial conditions..
x0 = 1
y0 = 1
y0p = 0
xn = 0.01
C
nout=11
open(nout,file='numout')
C
call Numerov(x0, xn, y0, y0p, numstp, funv, funu, xk, yk)
C
C ..write results..
do j=0,numstp,100
zk = yk(j)/( xk(j)**0.75 * exp(2/sqrt(xk(j))))
write(nout,*) xk(j), yk(j), zk
end do
C
stop 'All done'
end
C
C------------------------------------------------------------
C
function funu(x)
funu=0.
return
end
function funv(x)
funv=1./x**3
return
end
C
C------------------------------------------------------------
C
subroutine Numerov(x0, xn, y0, y0p, numstp, funv, funu, xk, yk)
C
dimension xk(0:numstp), yk(0:numstp)
C
C ..Numerov algorithm driver routine of the order of O(h**5)..
C for solving the second order differential equation
C
C d2y(x)/dx2=U(x)+V(x)*y(x)
C
C (see: Gonzales, Thompson, CIP, 11/1997)..
C
C INPUT:
C x0 -- starting value of x
C xn -- final value of x
C y0 -- starting value of y = y(x0)
C y0p -- starting value of dy/dx at x0
C y1 -- y1 = y(x1) is calculated only at the start
C of calculations
C funv - external routine funv(x)
C funu - external routine funu(x)
C
C OUTPUT:
C xk(0:numstp) -- independent variable x
C yk(0:numstp) -- dependent variable y=y(x)
C ..
external funv, funu
C
C ..calculate step size..
h = (xn-x0)/numstp
C
xk(0) = x0
yk(0) = y0
C
C ..then calculate y1..
y1 = y1Num(x0, y0, y0p, h, funv, funu)
xk(1) = x0 + h
yk(1) = y1
C
C ..and next all other steps..
do j=2,numstp
y2 = y2Num(x0, y0, y0p, h, y1, funv, funu, y2)
xk(j) = x0 + h
yk(j) = y2
x0 = x0 + h
y0 = y1
y1 = y2
end do ! end j
C
return
end
C
C------------------------------------------------------------
C
function y2Num(x0, y0, y0p, h, y1, funv, funu, y2)
C
C ..this subroutine gives the one step solution of the equation
C d2y(x)/dx2=U(x)+V(x)*y(x) using Numerov algorithm.
C The initial conditions are: y(x0) = y0, yprim(x0) = yp0
C The subroutine makes one step h in independent variable x
C giving the value y2 = y(x+2h) if the two previous values
C of y, y0 = y(x0) and y1 = y(x+h), are given
C (see: Gonzales, Thompson, CIP, 11/1997)..
C
C INPUT:
C x0 -- starting value of x
C y0 -- starting value of y = y(x0)
C y0p -- starting value of dy/dx at x0
C y1 -- y1 = y(x1) is calculated separately at the start
C of calculations
C
C OUTPUT:
C y2 - output value y2 = y(x2)
C
external funv, funu
C
x1 = x0 + h
x2 = x0 + 2*h
U0 = funu(x0)
U1 = funu(x1)
U2 = funu(x2)
V0 = funv(x0)
V1 = funv(x1)
V2 = funv(x2)
F0 = U0 + V0*y0
F1 = U1 + V1*y1
C
h2=h*h
y2Num = (2*y1 - y0 + h2/12*(U2 + 10*F1 + F0))/(1-V2*h2/12)
C
return
end
C
C------------------------------------------------------------
C
function y1Num(x0, y0, y0p, h, funv, funu)
C
C ..Calculate y1 = y(x1) for Numerov O(h**5) algorithm
C (see: Gonzales, Thompson, CIP, 11/1997)..
C
C INPUT:
C x0 -- starting value of x
C y0 -- starting value of y = y(x0)
C y0p -- starting value of dy/dx at x0
C
C OUTPUT:
C y1 = y(x+h)
C
external funv, funu
C
x1 = x0 + h
x2 = x1 + h
U0 = funu(x0)
U1 = funu(x1)
U2 = funu(x2)
V0 = funv(x0)
V1 = funv(x1)
V2 = funv(x2)
F0 = U0 + V0*y0
C
h2=h*h
a11 = 1 - V1*h2/4
a12 = V2*h2/24
a21 = -2 - 5*V1*h2
a22 = 1 - V2*h2/12
b1 = y0 + h*y0p + h2*(7*F0 + 6* U1 - U2)/24
b2 = -y0 + h2*(F0 + 10*U1+U2)/12
C
y1Num = (a22*b1 - a12*b2)/(a11 - a12*a21)
return
end
C
C------------------------------------------------------------
C