Fortran学习18：数值计算方法1

最新推荐文章于 2022-10-17 11:07:04 发布

鲁智深坐捻绣花针

最新推荐文章于 2022-10-17 11:07:04 发布

阅读量1.8k

点赞数 4

分类专栏： Fortran 文章标签：学习开发语言

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

版权

Fortran 专栏收录该内容

24 篇文章 24 订阅

订阅专栏

指针之后的问题，包括module、调用其他语言以及GUI，我后面有时间再做吧，感觉对于我来说没啥大的感觉，因为module他其实和c++的class是异曲同工的，没有多少需要用大力气学习的地方，至于调用其他语言，那就更算了，没见过谁Fortran调用c和VB的，感觉没啥用，要调用，也是python调用c++，或者后面看看Python能不能调用Fortran，至于GUI，算了算了，书上的gui画面难看死了，还不如抽时间去学一学Qt，或者是python的canvas、pyqt5和mayavi，再不济学学matlab作图Fortran的作图，就了解一下，知道他还能做gui就行了。

现在主要还是数值计算方法，这个非常重要，虽然Fortran他是又imsl库来解决的，可是你要是不知道数值计算如何进行的，那后面正反演方程你就不会写，更做不出来东西。

一、求解非线性函数：

主要是三种方法：1、二分法(Bisection)2、割线法(Secant)3、牛顿法(Newton,也是最重要的)

1、Bisection：

基本思想：

1、先猜两个值a,b，要让f(a)*f(b)<0，即根据零点定理，再[a,b]之间必有一点c使得f(c)=0

2、取 $c=\frac{a+b}{2}$ ,如果f(c)=0，就找到了正解，如果不等于零，计算f(c)，判断f(a)*f(c)以及f(c)*f(b)，看哪个小于零，在根据零点定理，重复1、2，直到f(c)=0或者是c的取值达到精度。

来一个 $f(x)=sin(x)-x$

real function func(input)
    implicit none
    real ::input
    func=input-sin(input)
    return
end function
program main
    implicit none
    REAL,PARAMETER::lim=0.0001
    REAL i,j,mid
    REAL,external::func
    ! [-1,1]
    i=1.0
    j=-1.0
    do while((func(i)==0).or.(i-j<lim))
        if(func(i)==0)then
        write(*,*)i
        end if
        exit
        if(func(j)==0)then
        write(*,*)j
        end if
        exit
        if(func(i)*func(j)<0)then
            mid=(i+j)/2
            if(func(mid)*func(i)<0)then
                i=mid
            else
                j=mid
            end if
        end if
    end do 
    write(*,*)mid
stop "over"
end program

这个程序的输出很简单，因为f(0)=0了，所以输出是这样的：

书上给了一个方程 $f(x)=(x+3)*(x-3)$ ，我们把函数func改一下进行计算：

real function func(input)
    implicit none
    real ::input
    ! func=input-sin(input)
    func=(input+3)*(input-3)
    return
end function

当然这是我的写法，之后我再来把书上给的标准二分法写一遍（这是个模板，改改就能用的那种）：

module numerical
    implicit none
    real,PARAMETER::zero=0.00001
    contains
        real function bisect(A,B,func)
        implicit none
        real A,B
        real mid
        real FA,FB,Fmid
        real,external::func
        mid=(A+B)/2
        Fmid=func(mid)
        do while(abs(Fmid)>zero)
            FA=func(A)
            FB=func(B)
            if(FA*FB<0)then
                B=mid
                mid=(A+B)/2.0
            else
                A=mid
                mid=(A+B)/2.0
            end if 
            Fmid=func(mid)
        end do
        bisect=mid
        return 
        end function 
        real function f1(X)
            implicit none
            real X
            f1=X-SIN(X)
        end function 
end module  
program main
    use numerical
    implicit none
    real A,B
    real ans
    do while(.true.)
        write(*,*)"input two number you guess"
        read(*,*)A,B
        if((f1(A)*f1(B)<0))exit 
        write(*,*)"wrong area"
    end do 
    ! A=-6.0
    ! B=5.0
    ans=bisect(A,B,f1)
    write(*,*)ans
    stop "over"
    end program

咋说呢，他这个程序有问题的，对于这个函数来说，输入的两个数如果是关于原点对称，那就会输出很快，如果输入的上下限不是关于原点对称的，那运行起来就很慢了。

计算起来就很慢。可见，这个二分法在求 $f(x)=0$ 的时候就不太好用喽。

割线法：

1、先选出两个猜测值a、b

2、画一条通过(a,f(a)),(b,f(b))的一条直线，令这条直线和x轴的交点为c，检查f(c)是否等于0，如果是，就找到了一个解，

3、f(c)!=0的时候，令b、c为两个猜测值a、b，重复上面的步骤。

        real function secant(a,b,f)
            implicit none
            real::a,b,c
            real,external::f 
            real::fa,fb,fc
            fa=f(a)
            fb=f(b)
            c=a-fa*(b-a)/(fb-fa)
            fc=f(c)
            do while(abs(fc)>zero)
                a=b
                b=c
                fa=f(a)
                fb=f(b)
                c=a-fa*(b-a)/(fb-fa)
                fc=f(c)
            end do 
            secant=c
        RETURN
        end function

把这段代码加入到上面的module numerical就可以再主程序中调用了，很方便

牛顿法：

1、先做一个猜想值a

2、以 $f'(a)$ 作为斜率，结果(a,f(a))作一条直线，令这条直线和x轴的交点为b，检查f(b)=0？

如果是就找到一个解，如果不是则3

3、f(b)！=0的时候，重新令b作为新的猜想值a，重复上述步骤。

咋说呢，主要就是导数求解问题，再就没啥可以说的了，这个做完了，我都可以用其他的语言来写一个newton，那最后就用python来一个吧(所有的Fortran代码结束)。

    real function newtow(a,f,df)
        implicit none
        real::a
        real,external::f
        real,external::df
        real::b,fb
        b=a-f(a)/df(a)
        fb=f(b)
        do while(abs(fb)>zero)
            a=b 
            b=a-f(a)/df(a)
            fb=f(b)
        end do 
        newtow=b
        return 
    end function
    real function dfunc(x)
        implicit none
        real::x
        dfunc=COS(x)
        return
    end function

针对的函数还是 $f(x)=sin(x)$ 把这两个函数加入到module numerical就ok，调用就行

来试一下结果：

ans=nweton(a,func,dfunc)

我没有改输入格式，所以还是输入了两个数，但是真正运算的就是a

这是完整的代码：

module numerical
    implicit none
    real,PARAMETER::zero=0.0001
    contains
        real function bisect(A,B,func)
        implicit none
        real A,B
        real mid
        real FA,FB,Fmid
        real,external::func
        mid=(A+B)/2
        Fmid=func(mid)
        do while(abs(Fmid)>zero)
            FA=func(A)
            FB=func(B)
            if(FA*FB<0)then
                B=mid
                mid=(A+B)/2.0
            else
                A=mid
                mid=(A+B)/2.0
            end if 
            Fmid=func(mid)
        end do
        bisect=mid
        return 
        end function 
        real function secant(a,b,f)
            implicit none
            real::a,b,c
            real,external::f 
            real::fa,fb,fc
            fa=f(a)
            fb=f(b)
            c=a-fa*(b-a)/(fb-fa)
            fc=f(c)
            do while(abs(fc)>zero)
                a=b
                b=c
                fa=f(a)
                fb=f(b)
                c=a-fa*(b-a)/(fb-fa)
                fc=f(c)
            end do 
            secant=c
        RETURN
        end function
        real function newtow(a,f,df)
        implicit none
        real::a
        real,external::f
        real,external::df
        real::b,fb
        b=a-f(a)/df(a)
        fb=f(b)
        do while(abs(fb)>zero)
            a=b 
            b=a-f(a)/df(a)
            fb=f(b)
        end do 
        newtow=b
        return 
    end function
    real function dfunc(x)
        implicit none
        real::x
        dfunc=COS(x)
        return
    end function
        real function  func(x)
            implicit none
            real::x
            func=SIN(x)
            return 
        end function
end module
    program main
        use numerical
        implicit none
        real::a,b
        real::ans
        write(*,*)"input two numbers"
        read(*,*)a,b
        ans=bisect(a,b,func)!二分法
        ans=secant(a,b,func)!割线法
        ans=newtow(a,func,dfunc)!牛顿法
        write(*,"('X=',F8.4)")ans
        stop "over"
    end

拿去用吧，记得叫我好心人哈哈，那么到这里就基本上结束了，明天要开始线性代数的使用和学习了。

最后，附上我的python代码和做的图，有点丑，别介意，嘿嘿。

import numpy as np
import matplotlib.pyplot as plt
import math

def df(x)->float:
    ans=math.cos(x)
    return ans
def func(x)->float:
    ans=math.sin(x)
    return ans
def newton(a)->float:
    b=a-func(a)/df(a)
    fb=func(b)
    zero=0.00001
    while abs(fb)>zero:
        a=b
        b=a-func(a)/df(a)
        fb=func(b)
    return b

arr=np.arange(-1,1,0.0001)
arr2=[]
for i in arr:
    arr2.append(newton(i))
fig=plt.figure()
plt.plot(arr,arr2)
plt.show()