3.3边值问题的直接积分
这节编写时用到了已知结果方程来得到第二个数值,方法不严谨
program main
implicit none
integer(8)::n,i,j,k
real(8),external::su,exact
real(8)::h,sp(1000),spi(1000),r(1000),haha(1000),m
open(unit=99,file="xable.csv")
open(unit=100,file="yable.csv")
open(unit=101,file="exaction.csv")
n=500
h=0.1
sp(1)=0
sp(2)=1-0.5*(h+2)*exp(-h)
sp(2)=0.995*sp(2)
do j=1,n+1,1
r(j)=(j-1)*h
write(99,*)r(j)
end do
do i=2,n,1
sp(i+1)=2*sp(i)-sp(i-1)&
+h*h/12*(su(r(i+1))+su(r(i-1))+10*su(r(i)))
end do
!!linear correction
m=(sp(n)-sp(n-100))/(r(n)-r(n-100))
do k=1,n+1,1
spi(k)=sp(k)-m*(k-1)*h
write(100,*)spi(k)
end do
do i=1,n+1,1
haha(i)=exact(r(i))
write(101,*)haha(i)
end do
end program main
!!ϕ(r) 满足的 Poisson 方程
function su(x)
implicit none
real(8)::x
real(8)::su
su=-0.5*x*dexp(-x)
return
end
function exact(y)
implicit none
real(8)::y
real(8)::exact
exact=1-0.5*(y+2)*dexp(-y)
return
end