FORTRAN+计算物理学学习日记（2）

利用Fortran编写数值微分函数，利用向前向后两点式以及五点式求微分。

以上节例题为模板，编写向前向后两点式如下：

!在例一的条件上，计算两点一次插值的微商，编写向前两点式,向后两点式
program main
    implicit none
    integer(4)::i,n
    real(4)::x1(7),y1(7),x(100),y(100),p(100),q(100),z(100),h
    x1(1)=0.0
    x1(2)=0.5
    x1(3)=1.0
    x1(4)=1.5
    x1(5)=2.0
    x1(6)=2.5
    x1(7)=3.0
    y1(1)=0.0
    y1(2)=0.4794
    y1(3)=0.8451
    y1(4)=0.9975
    y1(5)=0.9093
    y1(6)=0.5985
    y1(7)=0.1411
    n=6
    open(unit=101,file="xable.csv",status="new")
    open(unit=102,file="yable.csv",status="new")
    open(unit=103,file="pable.csv",status="new")
    open(unit=104,file="qable.csv",status="new")
    open(unit=105,file="hable.csv",status="new")
    h=3.0/100
    do i=1,100,1
        x(i)=h*i
        write(101,*) x(i)
        call interp6(x1,y1,x(i),n,y(i))
        write(102,*)y(i)
    end do
    do i=1,99,1
        p(i)=(y(i+1)-y(i))/h
        write(103,*) p(i)
    end do
    q(1)=(y(1)-0.00)/h
    write(104,*)q(1)   
    do i=2,100,1
        q(i)=(y(i)-y(i-1))/h
        write(104,*) q(i)
    end do   
    close(101)
    close(102)
    close(103)
    close(104)
    
    end program main
    
subroutine interpolyb(x1,x,n,yy)
integer k,j
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::x
real,intent(out)::yy(n+1)
yy(1)=1;
do k=2,n+1
    yy(1)=yy(1)*(x-x1(k))/(x1(1)-x1(k))
end do
yy(n+1)=1;
do k=1,n
    yy(n+1)=yy(n+1)*(x-x1(k))/(x1(n+1)-x1(k));
end do
do j=2,n
    yy(j)=1
    do k=1,j-1
        yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k))
    end do
    do k=j+1,n+1 
    yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k)); 
    end do
end do
return 
    end subroutine interpolyb
    
subroutine interp6(x1,y1,x,n,y)
real,intent(out)::y
real yy(n+1)
integer::i
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::y1(n+1)
real,intent(in)::x
y=0
call interpolyb(x1,x,n,yy)
do i=1,n+1
    y=y+yy(i)*y1(i)
end do
end subroutine interp6

    

绘制如下图片：

 

 五点式代码及绘制图像：

!在例一的条件上，改动部分函数，编写五点式方程
program main
    implicit none
    integer(4)::i,n
    real(4)::x1(7),y1(7),x(103),y(103),z(103),h
    x1(1)=0.0
    x1(2)=0.5
    x1(3)=1.0
    x1(4)=1.5
    x1(5)=2.0
    x1(6)=2.5
    x1(7)=3.0
    y1(1)=0.0
    y1(2)=0.4794
    y1(3)=0.8451
    y1(4)=0.9975
    y1(5)=0.9093
    y1(6)=0.5985
    y1(7)=0.1411
    n=6
    open(unit=101,file="xxable.csv",status="new")
    open(unit=102,file="yyable.csv",status="new")
    open(unit=103,file="hable.csv",status="new")
    h=3.0/100
    do i=1,103,1
        x(i)=h*i-2*h
        write(101,*) x(i)
        call interp6(x1,y1,x(i),n,y(i))
        write(102,*)y(i)
    end do
    close(101)
    close(102)
    !编写五点式
    do i=3,101,1
        z(i)=(y(i-2)-8*y(i-1)+8*y(i+1)-y(i+2))/(12*h)
        write(103,*)z(i)
    end do
    close(103)
    end program main
    
subroutine interpolyb(x1,x,n,yy)
integer k,j
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::x
real,intent(out)::yy(n+1)
yy(1)=1;
do k=2,n+1
    yy(1)=yy(1)*(x-x1(k))/(x1(1)-x1(k))
end do
yy(n+1)=1;
do k=1,n
    yy(n+1)=yy(n+1)*(x-x1(k))/(x1(n+1)-x1(k));
end do
do j=2,n
    yy(j)=1
    do k=1,j-1
        yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k))
    end do
    do k=j+1,n+1 
    yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k)); 
    end do
end do
return 
    end subroutine interpolyb
    
subroutine interp6(x1,y1,x,n,y)
real,intent(out)::y
real yy(n+1)
integer::i
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::y1(n+1)
real,intent(in)::x
y=0
call interpolyb(x1,x,n,yy)
do i=1,n+1
    y=y+yy(i)*y1(i)
end do
end subroutine interp6