A sample of toy FDM code

by Yi Zhang

Solving equation

u_t-u''=f(x,u)

f(x,u)=-2+2*(\sin(u)-\sin(x^2))

Time integration is full implicit, and during each time step Newton’s method is employed.

PROGRAM MAIN
!--------------------------------------------------
!Implicit time evolution of equation
!u_t-u''=f(x,u)
!f(x,u)=-2+r*(sin(u)-sin(x^2)),r=2
!--------------------------------------------------
! t=time points
! dt=time step
! h=neighboring nodes distance
! u=unknown
! x=node location
! n=No. of unknowns

!USE SUNPERF

IMPLICIT NONE
INTEGER(4) :: n
REAL(4) :: elapse,time(2),etime
REAL(8) :: t=0,dt=0.1
REAL(8) :: h
REAL(8),ALLOCATABLE :: u(:),x(:)
INTEGER(4) :: count,i

n=50
h=1/(DBLE(n)+1.)
ALLOCATE (x(n),u(n))
x=(/(i*h,i=1,n)/)
u=(/((i*h)**0.5,i=1,n)/)

OPEN(UNIT=100,FILE='data',STATUS='UNKNOWN',ACTION='&
     &WRITE')
DO i=1,n
   WRITE(100,*) t,x(i),u(i)
END DO

DO WHILE (t<1.0)    
CALL newton(n,h,dt,u,x)    
t=t+dt    
DO i=1,n       
WRITE(100,*) t,x(i),u(i)    
END DO END DO CLOSE(UNIT=100) 
elapse=etime(time) 
WRITE(*,*) 'elapsed:',elapse 
END PROGRAM MAIN 
!-------------------------------------------------- 
SUBROUTINE newton(n,h,dt,u,x)   
IMPLICIT NONE   
INTEGER(4),INTENT(IN) :: n   
REAL(8),INTENT(IN) :: x(n),h,dt   
REAL(8),INTENT(INOUT) :: u(n)   
REAL(8) :: r=2,c(n),error,b(n,1),b0(n,1),a(2,n)   
REAL(8) :: alpha=1.0,beta=-1.0   
INTEGER(4) :: info,i,iter   
c=u*h*h/dt   
iter=0   
error=1   
DO WHILE (error>1.e-13)
     iter=iter+1
     a(1,:)=-1
     a(2,:)=2+h*h/dt
     b(:,1)=(-2+r*(SIN(u)-SIN(x*x)))*h*h+c
     b(n,1)=b(n,1)+1
     CALL DSBMV('u',n,1,alpha,a,2,u,1,beta,b,1)
     b=-b
     a(2,:)=a(2,:)-r*COS(u)*h*h
     CALL DGTSV(n,1,a(1,2:),a(2,:),a(1,2:),b,n,info)
     error=MAXVAL(ABS(b))
     u=u+b(:,1)
  END DO
  WRITE(*,*) 'INFO',info,'Number of iterations:',iter
END SUBROUTINE newton
!--------------------------------------------------

Output would be like:

% f95 -dalign implicit_newton.f95 -xlic_lib=sunperf
% a.out
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 4
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 INFO 0 Number of iterations: 3
 elapsed: 0.061407

Of course the exact solution is u(x)=x^2. Here is time evolution plot of the numerical solution.

implicit_newton

Advertisements