Shooting P1, Fortran

Shooting P1

Program Finite_Shooting
    Implicit None
     Real::aa=0.0,bb=atan(1.0),h,w(0:15),p,q,r,x,ex,&
          a(15),b(15),d(15),l(15),u(15),z(15),c(15)
    Real::u1(100,0:100),v(100,0:100),w1(100,0:100),k(4,2),g(4,2)
    integer::n=13,i

h=(bb-aa)/(n+1)

CALL SM(x,w1)

Write(*,*)"         x             exact              shooting      "
    Write(*,*)"_______________________________"
     Do i=0,n+1
        x=aa+i*h
        Write(*,22)x,ex(x),w1(1,i)
    end do

22 FORMAT(4(2X,f13.6))
End Program

Subroutine SM(x,w1)
Implicit None
    Real::aa=0.0,bb=atan(1.0),u(100,0:100),v(100,0:100),&
          w1(100,0:100),a,b,x,k(4,2),g(4,2),h,p,q,r,ex
 Integer::n=14,i
 a=1.0
 b=2.0
 h=(bb-aa)/n
 u(1,0)=a
 u(2,0)=0
 v(1,0)=0
 v(2,0)=1

Do i=0,n-1
    x=a+i*h
    k(1,1)=h*u(2,i)
    k(1,2)=h*(p(x)*u(2,i)+q(x)*u(1,i)+r(x))
    k(2,1)=h*(u(2,i)+0.5*k(1,2))
    k(2,2)=h*(p(x+(h/2.))*(u(2,i)+0.5*k(1,2))+q(x+(h/2.))*(u(1,i)+0.5*k(1,1))+&
    r(x+(h/2.)))
    k(3,1)=h*(u(2,i)+0.5*k(2,2))
    k(3,2)=h*(p(x+(h/2.))*(u(2,i)+0.5*k(2,2))+q(x+(h/2.))*(u(1,i)+0.5*k(2,1))+&
    r(x+(h/2.)))
    k(4,1)=h*(u(2,i)+k(3,2))
    k(4,2)=h*(p(x+h)*(u(2,i)+k(3,2))+q(x+h)*(u(1,i)+k(3,1))+r(x+h))
    u(1,i+1)=u(1,i)+(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1))/6.0
    u(2,i+1)=u(2,i)+(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2))/6.0

g(1,1)=h*v(2,i)
    g(1,2)=h*(p(x)*v(2,i)+q(x)*v(1,i))
    g(2,1)=h*(v(2,i)+0.5*g(1,2))
    g(2,2)=h*(p(x+(h/2.))*(v(2,i)+0.5*g(1,2))+q(x+(h/2.))*(v(1,i)+0.5*g(1,1)))
    g(3,1)=h*(v(2,i)+0.5*g(2,2))
    g(3,2)=h*(p(x+(h/2.))*(v(2,i)+0.5*g(2,2))+q(x+(h/2.))*(v(1,i)+0.5*g(2,1)))
    g(4,1)=h*(v(2,i)+g(3,2))
    g(4,2)=h*(p(x+h)*(v(2,i)+g(3,2))+q(x+h)*(v(1,i)+g(3,1)))
    v(1,i+1)=v(1,i)+(g(1,1)+2*g(2,1)+2*g(3,1)+g(4,1))/6.0
    v(2,i+1)=v(2,i)+(g(1,2)+2*g(2,2)+2*g(3,2)+g(4,2))/6.0
end do

do i=0,n
    w1(1,i)=u(1,i)+((b-u(1,n))*v(1,i))/v(1,n)
end do

End Subroutine

Real function p(x)
    implicit none
    real::x
    p=-2.0/x
    return
    end function p

Real function q(x)
    implicit none
    real::x
    q=2.0/(x**2.0)
    return
    end function q

Real function r(x)
    implicit none
    real::x
    r=sin(log(x))/(x**2.0)
    return
    end function r

Real function ex(x)
    implicit none
    real::x
    ex=-0.03920701320*x+(1.1392070132)*(1/x**2)-(3/10)*sin(log(x))-(1/10)*cos(log(x))
    return
    end function ex


run \| edit \| history \| help	0