谁的借势营销做的最好

问题:冈本:今益求精,头头事到\nNokia:赢得优良\nMotorola:恭喜中国足球,今天看自身的\nHUAWEI科技(science and technology):只要有1%的或许,大家就有99%的信教\nAlienware:追寻必须要经过的路,就是承受14亿梦想,毫不退让\n安踏体育:有时翻脸比不上翻身\n恰恰食物:国“仁”骄傲\n……\n

求解难题

Write a program to solve the following ordinary differential
equation by

  • basic Euler method
  • improved Euler method
  • four-order Runge-Kutta method

图片 1=3%20%5Cend%7Bcases%7D%20x%5Cin%5B0,%201.5%5D)

and calculate y(1.5) with stepsize=0.1, 0.1/2, 0.1/4, 0.1/8

Compare it with analytic solution (in figure)

图片 2%20=%20%5Cdfrac%7B3%7D%7B1+x%5E3%7D)

回答:

主程序

program main
    use ODE
    implicit none
    real :: x0=0.0,y0=3.0,a=0.0,b=1.5,h=0.1
    integer :: n,i,j
    character(len=512) :: filename


    print *,'The analytic solution of y(1.5)= ',analytic_solution(1.5)
    print *

    print *,'Basic Euler Method'
    print '(5x,a,f3.1,3x,a,f3.1,3x,a,f3.1,3x,a,f3.1)','Set initial  x0= ',x0,'y0= ',y0,'a= ',a,'b= ',b
    do j=0,3
        h=0.1/2**j
        print '(5x,a,f6.4)','Set h= ',h
        call init(x0,y0,a,b,h)
        n=size(x)-1
        call Basic_Euler
        print '(8x,a,f10.8)','y(1.5)= ',y(n)
        write(filename, *) j
        filename='BEM_'//Trim(AdjustL(filename))//'.txt'
        open(101,file=filename)
        write(101,'(f3.1,f10.8)') (x(i),y(i),i=0,n)
        close(101)
        deallocate(x,y)
    end do
    print *

    print *,'Improved Euler Method'
    print '(5x,a,f3.1,3x,a,f3.1,3x,a,f3.1,3x,a,f3.1)','Set initial  x0= ',x0,'y0= ',y0,'a= ',a,'b= ',b
    do j=0,3
        h=0.1/2**j
        print '(5x,a,f6.4)','Set h= ',h
        call init(x0,y0,a,b,h)
        n=size(x)-1
        call Improved_Euler
        print '(8x,a,f10.8)','y(1.5)= ',y(n)
        write(filename, *) j
        filename='IEM_'//Trim(AdjustL(filename))//'.txt'
        open(101,file=filename)
        write(101,'(f3.1,f10.8)') (x(i),y(i),i=0,n)
        close(101)
        deallocate(x,y)
    end do
    print *

    print *,'Four Order Runge-Kutta Method'
    print '(5x,a,f3.1,3x,a,f3.1,3x,a,f3.1,3x,a,f3.1)','Set initial  x0= ',x0,'y0= ',y0,'a= ',a,'b= ',b
    do j=0,3
        h=0.1/2**j
        print '(5x,a,f6.4)','Set h= ',h
        call init(x0,y0,a,b,h)
        n=size(x)-1
        call Four_Order_Runge_Kutta
        print '(8x,a,f10.8)','y(1.5)= ',y(n)
        write(filename, *) j
        filename='RKM_'//Trim(AdjustL(filename))//'.txt'
        open(101,file=filename)
        write(101,'(f3.1,f10.8)') (x(i),y(i),i=0,n)
        close(101)
        deallocate(x,y)
    end do
    print *

end program main

求解ODE方程的要害办法写在ODE模块中:

module ODE
    implicit none
    private
    public :: x,y,init,Basic_Euler,Improved_Euler,Four_Order_Runge_Kutta,analytic_solution

    real,allocatable :: x(:),y(:)
    real :: x0,y0,a,b,h
    integer :: n,i

contains

    function f(x,y)
        implicit none
        real :: f,x,y
        f=-x*x*y*y
    end function f

    function analytic_solution(x) result(f)
        implicit none
        real :: f,x
        f=3.0/(1+x*x*x)
    end function

    subroutine init(x0_,y0_,a_,b_,h_)
        implicit none
        real :: x0_,y0_,a_,b_,h_
        x0=x0_
        y0=y0_
        a=a_
        b=b_
        h=h_
        n=int((b-a)/h)
        allocate(x(0:n),y(0:n))
        x=(/ (a+i*h,i=0,n) /)
        y=0
        y(0)=y0
    end subroutine init

    subroutine Basic_Euler()
        implicit none
        do i=1,n
            y(i)=y(i-1)+h*f(x(i-1),y(i-1))
        end do
    end subroutine Basic_Euler

    subroutine Improved_Euler()
        implicit none
        real :: y_
        do i=1,n
            y_=y(i-1)+h*f(x(i-1),y(i-1))
            y(i)=y(i-1)+h/2*(f(x(i-1),y(i-1))+f(x(i),y_))
        end do
    end subroutine Improved_Euler

    subroutine Four_Order_Runge_Kutta()
        implicit none
        real :: k1,k2,k3,k4
        do i=1,n
            k1=f(x(i-1),y(i-1))
            k2=f(x(i-1)+h/2,y(i-1)+h/2*k1)
            k3=f(x(i-1)+h/2,y(i-1)+h/2*k2)
            k4=f(x(i-1)+h,y(i-1)+h*k3)
            y(i)=y(i-1)+h/6*(k1+2*k2+2*k3+k4)
        end do
    end subroutine Four_Order_Runge_Kutta
end module ODE

其中init方法是用来初叶化:

    subroutine init(x0_,y0_,a_,b_,h_)
        implicit none
        real :: x0_,y0_,a_,b_,h_
        x0=x0_
        y0=y0_
        a=a_
        b=b_
        h=h_
        n=int((b-a)/h)
        allocate(x(0:n),y(0:n))
        x=(/ (a+i*h,i=0,n) /)
        y=0
        y(0)=y0
    end subroutine init

analytic_solution是数值解:

    function analytic_solution(x) result(f)
        implicit none
        real :: f,x
        f=3.0/(1+x*x*x)
    end function

不是广告,小编只是感觉本次大家合营社的文案不错
图片 3

Basic Euler Method

流程图:

图片 4

原理:
![][3]
[3]:
http://latex.codecogs.com/gif.latex?\begin{cases}%20y(x_{n+1})=y(x_n)+hf(x_n,y(x_n))+O(h^2)%20\%20y(x_0)=y_0%20\end{cases}

代码:

    subroutine Basic_Euler()
        implicit none
        do i=1,n
            y(i)=y(i-1)+h*f(x(i-1),y(i-1))
        end do
    end subroutine Basic_Euler

出口结果:

图片 5

Basic Euler 方法

图例:

图片 6

区别h每步迭代结果

Improved Euler Method

流程图:

图片 7

原理:
![][4]
[4]:
http://latex.codecogs.com/gif.latex?\begin{cases}%20\bar%20y_{n+1}=y_n+hf(x_n,y_n)%20\%20y(x_{n+1})=y(x_n)+\dfrac{h}{2}[f(x_n,y(x_n))+f(x_{n+1},\bar%20y_{n+1})]+O(h^2)%20\%20y(x_0)=y_0%20\end{cases}

代码:

    subroutine Improved_Euler()
        implicit none
        real :: y_
        do i=1,n
            y_=y(i-1)+h*f(x(i-1),y(i-1))
            y(i)=y(i-1)+h/2*(f(x(i-1),y(i-1))+f(x(i),y_))
        end do
    end subroutine Improved_Euler

出口结果:

图片 8

Improved Euler 方法

图例:

图片 9

差别h每步迭代结果

Four Order Runge-Kutta Method

流程图:

图片 10

原理:
![][5]
[5]:
http://latex.codecogs.com/gif.latex?\begin{cases}%20y_{n+1}=y_n+\dfrac{h}{6}(K_1+2K_2+2K_3+K_4)%20\%20K_1=f(x_n,y_n)%20\%20K_2=f(x_n+\dfrac{h}{2},y_n+\dfrac{h}{2}K_1)%20\%20K_3=f(x_n+\dfrac{h}{2},y_n+\dfrac{h}{2}K_2)%20\%20K_4=f(x_n+h,y_n+hK_3)%20\end{cases}

代码:

    subroutine Four_Order_Runge_Kutta()
        implicit none
        real :: k1,k2,k3,k4
        do i=1,n
            k1=f(x(i-1),y(i-1))
            k2=f(x(i-1)+h/2,y(i-1)+h/2*k1)
            k3=f(x(i-1)+h/2,y(i-1)+h/2*k2)
            k4=f(x(i-1)+h,y(i-1)+h*k3)
            y(i)=y(i-1)+h/6*(k1+2*k2+2*k3+k4)
        end do
    end subroutine Four_Order_Runge_Kutta

输出结果:

图片 11

4 Order Runge-Kutta 方法

图例:

图片 12

不一样h每步迭代结果

二种艺术结果比较

图片 13

h=0.1时

图片 14

h=0.1/8时

能够看到,当h非常大时,三种艺术的分裂依旧比非常的大的,当h慢慢减时辰,两种方法的结果已基本同样。

网站地图xml地图