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.

### 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.