Fortran学习21：数值计算方法4

最新推荐文章于 2024-09-26 06:59:54 发布

鲁智深坐捻绣花针

最新推荐文章于 2024-09-26 06:59:54 发布

阅读量1k

点赞数 4

分类专栏： Fortran 文章标签：学习算法

本文链接：https://blog.csdn.net/faltas/article/details/127349008

版权

Fortran 专栏收录该内容

24 篇文章 26 订阅

订阅专栏

今天要学习的就是插值方法了，这一节量很大的，我感觉今天一天搞不完的，能学到那就算哪吧。

数值计算方法我建议可以取读一本书，高尚、秦斌合写的计算方法，算是对母校的书籍进行一个推广吧。

1、拉格朗日插值方法：

经过n个数据点，生成最高项为n-1的多项式，

假设：又(X1,Y1)、(X2,Y2)、(X3,Y3)......(XN,YN)个数据连，经过拉格朗日插值法获得的多项式就是：

$g(x)=\frac{(x-x_{2})(x-x_{3})(x-x_{4})...(x-x_{n})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})...(x_{1}-x_{n})}y_{1}+ \frac{(x-x_{1})(x-x_{3})(x-x_{4})...(x-x_{n})}{(x_{2}-x_{1})(x_{2}-x_{3})(x_{2}-x_{4})...(x_{2}-x_{n})}y_{2}+....+ \frac{(x-x_{1})(x-x_{2})(x-x_{3})...(x-x_{n})}{(x_{n}-x_{1})(x_{n}-x_{2})(x_{n}-x_{3})...(x_{n}-x_{n-1})}y_{n}$

我无语的地方是这个用的是作者自己搞的sgl库来绘制图形，emm...咋说呢，我没有光驱，那就算了，直接写代码吧，后面用matlab来作图。

module interpolate_utility
    implicit none
    type point
        real x,y
    end type 
    real,parameter::pi=3.1415926
    real,PARAMETER::xmin=0.0,xmax=pi*3.0
    integer,PARAMETER::n=10,np=100
    type(point)::datas(n)
    type(point)::interpolate(np)
    contains
    subroutine generateData(func)
        implicit none
        real,external::func
        real r,width
        integer i
        width=(xmax-xmin)/(n-1)
        r=0
        do i=1,n
            datas(i)%x=r
            datas(i)%y=func(r)
            r=r+width
        end do
    end subroutine
    real function lagrange(x)
        implicit none
        real x
        real coeff
        integer i,j
        lagrange=0
        do i=1,n
            coeff=1
            do j=1,n
            if(i/=j)coeff=coeff*(x-datas(j)%x)/(datas(i)%x-datas(j)%x)
            end do 
        lagrange=lagrange+coeff*datas(i)%y
        end do 
    end function
end module
program main
    use interpolate_utility
    implicit none
    REAL,intrinsic:: sin
    real xinc,x 
    integer i
    integer,PARAMETER:: outinterface=10
    CHARACTER(len=20)::filename="plot.txt"
    CALL generateData(sin)
    x=0
    xinc=(xmax-xmin)/(np-1)
    do i=1,np
        interpolate(i)%x=x
        interpolate(i)%y=lagrange(x)
        x=x+xinc
    end do
    open(unit=outinterface, file=filename)
    do i=1,np
        write(outinterface,*)interpolate(i)%x,interpolate(i)%y
    end do 
    
    close(outinterface)
    stop "over"
    end program

Fortran计算起来速度还是快，用来写原理了，至于图像，我觉得python或者matlab都比Fortran要好得多呢，所以我试试把他们写在一个txt文件里面，然后用python读出来，做个图。

import pandas as pd
import matplotlib.pyplot as plt
filename="H:/研究生相关文献和程序，代码等/code/fortranproject/plotla.txt"
data = pd.read_table(filename,header=None,delim_whitespace=True)
x=data[0]
y=data[1]
fig=plt.figure()
plt.plot(x,y)
plt.show()

啊哈哈，真漂亮。

2、牛顿法：

lagrange有个最大的缺点，每一次要插入新数值的时候，就要把全部的数据点代入，重新进行计算，很麻烦，可以用牛顿法，牛顿法事先把部分计算创建成表格，，插值的时候引用这个表格里的数值来加快计算。

$g(x)=\sum_{n=0}^{k}\binom{s}{n}\Delta ^{n}f_{0}$

其中：

\binom{s}{n}

$\binom{s}{0}=1$

$\binom{s}{1}=s$

$\binom{s}{2}=\frac{s(s-1)}{2!}$

$\binom{s}{2}=\frac{s(s-1)(s-2)(s-3)...(s-n+1)}{n!}$

s的定义是：

$s=\frac{x-x_{0}}{h}$

h是数据点之间的距离，使用这个方法的时候，数据点之间的距离必须相等，即

$x_{1}-x_{0}=x_{2}-x_{1}=h$

$\Delta ^{n}f_{i}$ :

$\Delta ^{n}f_{i}=\Delta ^{n-1}f_{i+1}-\Delta ^{n-1}f_{i}$

具体的表格是这样子的：

I	$f_{i}$	$\Delta ^{1}f_{i}$	$\Delta ^{2}f_{i}$	$\Delta ^{3}f_{i}$
1
2
3
4

    subroutine buildtable()
        implicit none
        integer row,col,i
        table=0
        do i=1,n
            table(i,1)=datas(i)%y
        end do
        do col=2,n
            do row=1,n-col+1
                table(col,row)=table(row+1,col-1)-table(row,col-1)
            end do
        end do
    end subroutine 

    real function newton(x,th,num)
        implicit none
        real x
        integer th,num
        real s,sum,coeff
        integer i,j,f
        if(th+num-1>n)then
            write(*,*)"data is empty"
        return
        end if
        newton=table(th,1)
        s=(x-datas(th)%x)/width
        f=1
        coeff=1.0
        do i=1,num-1
            f=f*i
            coeff=coeff*(s-i+1)
            newton=newton+coeff*table(th,i+1)/f
        end do
        return
    end function

调用就行了。