Fortran学习20：数值计算方法3

今天来看看微积分，什么是微积分，好多人会给出自己的理解，我的理解就是，一个非规则几何体、面，太难算了，就把他给分，分成足够多的小块，然后就发现，每一个小块的面积相加就是整体的面积，分的过程就是微分，相加的过程就是积分。

用数值方法来积分的时候，主要是应用在函数不存在的时候，我们从实验中得到的一连串数据，没有办法刚好使用一些已知的函数来表示，那么这个时候就用其他的算法来计算了。

1、梯形法求积分：

这个是最简单的，使用梯形法，把所需要积分的图形，用许多的小梯形方块填满，这样就近似求解了。

module integeral
    implicit none
    real,PARAMETER::PI=3.1415926
    contains
        subroutine GenerateData(datas,width,func)
            ! 生成数列
            real datas(:),width
            real,external::func
            real r
            integer i,n
            n=size(datas, dim=1)
            width=PI/(n-1)
            r=0
            do i=1,n
                datas(i)=func(r)
                r=r+width
            end do     
        end subroutine
        real function Trape_Integral(datas,width)
            implicit none
            real datas(:)
            real::width,sum
            ! sum：计算所有上底加下底除以2的和
            integer i,n
            n=size(datas, dim=1)
            sum=(datas(1)+datas(n))/2.0
            ！（上底加下底）/2
            do i=2,n-1
                sum=sum+datas(i)
                ! 累加边长
            end do 
            Trape_Integral=sum*width
            ! 计算面积和
            return 
        end function
end module
program main
    use integeral
    implicit none
    integer,PARAMETER::n=10
    real::datas(n),width
    real ans
    real,intrinsic::sin
    ! 仿真用来生成数据的函数
    CALL GenerateData(datas,width,sin)
    ans= Trape_Integral(datas,width)
    write(*,"('ans=',F4.2)")ans
    stop "over"
    end program

program main中的n是精度，n越大，拟合的效果就越好 

上面是n=10的情况，再看看n=10000：

 n在这里就是微分的精度，微分的精度越高，拟合的效果就越好

辛普森法求积分：

辛普森法的原理就是把函数图像用很多段的二次曲线来近似图形，在计算这些二次曲线所形成的面积，但是使用辛普森法有个条件：必须要是奇数个数据才能够计算，因为每一条二次曲线都要从数据中取三个点，每一条二次曲线的积分公式：

把每一段二次曲线累加起来，就可以得到整个辛普森法积分的公式：

上代码：

        real function Simpson(datas,width)
            implicit none
            real datas(:),width
            real sum
            integer i,n
            n=size(datas, dim=1)
            if(mod(n,2)==0)then
                write(*,*)"you must enter an odd number"
            stop "error number"
            end if
            sum=datas(1)+datas(n)
            ! 先算出头尾的和
            do i=2,n-1
                if(mod(i,2)==0)then
                    sum=sum+4*datas(i)
                else
                    sum=sum+2*datas(i)
                end if
            end do 
            Simpson=sum*width/3.0
            RETURN
        end function 

把这一段加到上面的module里面，直接调用就行，就是在调用的时候，program main中的n必须是奇数，否则会在第七行报错。

上述两种方法算的数值一一模一样的，就不贴代码了。

这里介绍三种Fortran计算时间的办法：

1、cpu_time()

real::time_start,time_end,time
call cpu_time(time_start)
......
call cpu_time(time_end)
time=time_start-time_end

这个只用于单线程，而且数值必须要是real，如果使用在多线程中，Fortran并不会报错，他输出的实际是所有线程的时间和

2、DATE_AND_TIME(dat1,dat2)

CHARACTER(len=20)::dat1,dat2
CALL DATE_AND_TIME(dat1,dat2)

输出的dat1就是当天的日期

3、system_clock(time)

这个输出的是毫秒，多线程用的是这个方法

integer time
CALL system_clock(time)

好了，今天就到这里吧，要去看论文了。。。