3.5打靶法
边值问题
program main
implicit none
integer(8)::n,i
real(8)::h,error,delta,ddelta,xs(100),yr(100),yr1(100),ka(100),kb(100),yn(100),ynn(100),re(100),pi
real(8),external:: g5
pi=3.1415926
n=40
h=1.0/n
error=10e-5
delta=0.8
ddelta=0.1
open(unit=100,file="xable.csv")
open(unit=101,file="numerov.csv")
open(unit=102,file="exact.csv")
do i=1,n+1,1
xs(i)=h*(i-1)
write(100,*)xs(i)
end do
close(100)
!! The Numerov Algorithm
call interp(delta,h,n,ynn)
do while(ddelta>error)
delta=delta-ddelta
call interp(delta,h,n,yn)
if((yn(n+1)-1)*(ynn(n+1)-1)>0)then
else
delta=delta+ddelta
ddelta=ddelta/2.0
end if
end do
do i=1,n+1,1
write(101,*)yn(i)
end do
close(101)
!!The exact result
do i=1,n+1,1
re(i)=g5(xs(i))
write(102,*)re(i)
end do
close(102)
end program
function g5(x)
implicit none
real(8)::x,pi
real(8)::g5
pi=3.1415926
g5=cos(1.0/2*pi*x)+2.0*sin(1.0/2*pi*x)-1.0
return
end
subroutine interp(delta,h,n,yx)
integer(8)::i
real(8)::con,pi
real(8),intent(in)::delta
real(8),intent(in)::h
integer(8),intent(in)::n
real(8),intent(out)::yx(100)
yx(1)=0
yx(2)=delta
pi=3.1415926
con=h*h*pi*pi/48.0
do i=2,n,1
yx(i+1)=2.0*(1.0-5*con)*yx(i)-(1+con)*yx(i-1)-h*h*pi*pi/4.0
yx(i+1)=yx(i+1)/(1+con)
end do
return
end subroutine interp
本征值问题
program main
implicit none
integer(8)::i,n
real(8)::h,delta,tol,k,dk,x(500),y(500),p1(500),p2(500),a,pp(500),p(100)
real(8),external::g4
n=50
h=1.0/n
delta=0.6
tol=10**(-5)
k=12.0
dk=1.0
open(unit=100,file="xable.csv")
open(unit=101,file="normalization.csv")
open(unit=102,file="exact.csv")
call interp(k,h,n,delta,p1)
do while(abs(dk)>tol)
k=k+dk
call interp(k,h,n,delta,p2)
if(p1(n+1)*p2(n+1)>0)then
else
k=k-dk
dk=dk/2.0
end if
end do
do i=1,n+1,1
x(i)=(i-1)*h
write(100,*)x(i)
end do
!!the normalization
a=0
do i=1,n,1
a=a+h/2*(p2(i)*p2(i)+p2(i+1)*p2(i+1))
end do
do i=1,n+1,1
p2(i)=p2(i)/dsqrt(a)
write(101,*)p2(i)
end do
!!exact solution
do i=1,n+1,1
y(i)=g4(x(i))
write(102,*)y(i)
end do
close(100)
close(101)
close(102)
end program main
function g4(x)
implicit none
real(8)::x,pi
real(8)::g4
pi=3.1415926
g4=sqrt(2.0)*dsin(4.0*pi*x)
return
end
subroutine interp(k,h,n,delta,p)
implicit none
real(8),intent(in)::k,h,delta
integer(8),intent(in)::n
real(8),intent(out)::p(100)
real(8) con
integer(8) j
p(1)=0
p(2)=delta
con=(k*h)*(k*h)/12.0
do j=2,n,1
p(j+1)=2.0*(1-5*con)*p(j)-(1+con)*p(j-1)
p(j+1)=p(j+1)/(1+con)
end do
end subroutine interp