Mesh Mess

FEM, CFD, CAE, and all that.

Tag: Finite difference

What’s the difference?

Today I was asked for a one-sentence comment of the difference between finite difference method and finite element method. What I came up was:

In finite difference method we approximate the differential expression, or the partial derivative, while in finite element method we approximate the space in which the solution lies.

Hopefully this shed some light on the issue.

Advertisements

A sample of toy FDM code

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