implicitscheme.f90

求解粘性Burger方程和非粘性Burger方程的各种差分格式

subroutine SchemeTwo(uOne,uZero,lenX,deltX,deltT,coeff)
! The difference scheme on page 140
! Lag Nonlinear Term
! coeff is with the second order term

!	integer,parameter :: maxLen=20
	real uOne(0:lenX),uZero(0:lenX)
	integer k
	real cofFir,cofSec          ! first order and second order's coeff in the diffition eqution 
	real a(2:lenX-1),b(1:lenX-1),c(1:lenX-2)

	cofSec=coeff*deltT/deltX**2;    b(1)=2*cofSec+1
	cofFir=deltT*uZero(1)/(2*deltX);c(1)=cofFir-cofSec

	do k=2,lenX-2
		cofFir=deltT*uZero(k)/(2*deltX)
		a(k)=-(cofFir+cofSec)
		b(k)=2*cofSec+1
		c(k)=cofFir-cofSec
	end do

	cofFir=deltT*uZero(lenX-1)/(2*deltX);a(lenX-1)=-(cofFir+cofSec)
	b(lenX-1)=2*cofSec+1

	call ThomasAlgorithm(a,b,c,lenX-1,uZero(1:lenX-1))

	uOne(1:lenX-1)=uZero(1:lenX-1)

	uOne(0)=bounLeft(time)                               ! Dirichlet boundary value
	uOne(lenX)=bounRight(time)

end subroutine SchemeTwo

!================================================================================================
subroutine SchemeThree(uOne,uZero,lenX,deltX,deltT,coeff)
! The difference scheme on page 140
! Linearize About previous Time Step
! coeff is with the second order term

!	integer,parameter :: maxLen=20
	real uOne(0:lenX),uZero(0:lenX)
	integer k
	real cofFir,cofSec          ! first order and second order's coeff in the diffition eqution 
	real a(2:lenX-1),b(1:lenX-1),c(1:lenX-2),r(1:lenX-1)

	cofSec=coeff*deltT/deltX**2  
	cofFir=deltT/deltX;

	b(1)=1+2*cofSec+cofFir/2.*(uZero(2)-uZero(0))
	c(1)=cofFir/2*uZero(1)-cofSec

	do k=2,lenX-2
		a(k)=-cofFir/2.*uZero(k)-cofSec
		b(k)=1+2*cofSec+cofFir/2.*(uZero(k+1)-uZero(k-1))
		c(k)=cofFir/2*uZero(k)-cofSec
	end do

	a(lenX-1)=-cofFir/2.*uZero(lenX-1)-cofSec
	b(lenX-1)=1+2*cofSec+cofFir/2.*(uZero(lenX)-uZero(lenX-2))

	do k=1,lenX-1
		r(k)=-cofFir/2.*(uZero(k+1)-uZero(k-1))+cofSec*(uZero(k+1)-2*uZero(k)+uZero(k-1))
	end do

	call ThomasAlgorithm(a,b,c,lenX-1,r)

	uOne(1:lenX-1)=r+uZero(1:lenX-1) 

	uOne(0)=bounLeft(time)                               ! Dirichlet boundary value
	uOne(lenX)=bounRight(time)

end subroutine SchemeThree

!===============================================================================================
subroutine SchemeSix(uOne,uZero,lenX,deltX,deltT,coeff,time)
! The difference scheme on page 141
! Newton's Method
! coeff is with the second order term
	use comNumFun
	real uOne(0:lenX),uZero(0:lenX)
	integer k
	real :: cofFir,cofSec,err         ! first order and second order's coeff in the diffition eqution 
	real a(2:lenX-1),b(1:lenX-1),c(1:lenX-2),r(1:lenX-1)
	
	err=1
	cofSec=coeff*deltT/deltX**2  
	cofFir=deltT/deltX
    uOne(0:lenX)=uZero(0:lenX)
    do while( abs(err)>1e-6 ) 
		b(1)=1+2*cofSec+cofFir/2.*(uOne(2)-uOne(0))
		c(1)=cofFir/2.*uOne(1)-cofSec

		do k=2,lenX-2
			a(k)=-cofFir/2.*uOne(k)-cofSec
			b(k)=1+2*cofSec+cofFir/2.*(uOne(k+1)-uOne(k-1))
			c(k)=cofFir/2*uOne(k)-cofSec
		end do

		a(lenX-1)=-cofFir/2.*uOne(lenX-1)-cofSec
		b(lenX-1)=1+2*cofSec+cofFir/2.*(uOne(lenX)-uOne(lenX-2))

		do k=1,lenX-1
			r(k)=-uOne(k)-cofFir/2.*(uOne(k+1)-uOne(k-1))+cofSec*(uOne(k+1)-2*uOne(k)+uOne(k-1))+uZero(k)
		end do

		call ThomasAlgorithm(a,b,c,lenX-1,r)
		err=maxArr(r,lenX-1)	
		uOne(1:lenX-1)=r+uOne(1:lenX-1)
		uOne(0)=bounLeft(time)                               ! Dirichlet boundary value
		uOne(lenX)=bounRight(time)		
	end do 

end subroutine SchemeSix

!================================================================================================

subroutine SchemeSeven(uOne,uZero,lenX,deltX,deltT,coeff,time)
! The difference scheme on page 347 (6.10.7)
! Newton's Method
! coeff is with the second order term
	use comNumFun
	real uOne(0:lenX),uZero(0:lenX)
	integer k
	real :: cofFir,cofSec,err         ! first order and second order's coeff in the diffition eqution 
	real a(2:lenX-1),b(1:lenX-1),c(1:lenX-2),r(1:lenX-1)
	real deltU(1:lenX-1)
		
	err=1
	cofSec=coeff*deltT/deltX**2  
	cofFir=deltT/deltX
    uOne(0:lenX)=uZero(0:lenX)
    do while( abs(err)>1e-3 )
		do k=1,lenX-1
			if(uOne(k)<0) then
				deltU(k)=(uOne(k+1)-uOne(k))/deltX
			else
				deltU(k)=(uOne(k)-uOne(k-1))/deltX
			end if
		end do

		if(uOne(1)<0) then 
			b(1)=1+2*cofSec+deltU(1)*deltT-cofFir*uOne(1)
			c(1)=cofFir*uOne(1)-cofSec
		else 
			b(1)=1+2*cofSec+deltU(1)*deltT+cofFir*uOne(1)
			c(1)=-cofSec
		end if			

		do k=2,lenX-2
			if (uOne(k)<0) then
				a(k)=-cofSec
				b(k)=1+2*cofSec+deltU(k)*deltT-cofFir*uOne(k)
				c(k)=cofFir*uOne(k)-cofSec
			else 
				a(k)=-cofSec-cofFir*uOne(k)
				b(k)=1+2*cofSec+deltU(k)*deltT+cofFir*uOne(k)
				c(k)=-cofSec
			end if
		end do

		if(uOne(lenX-1)<0) then
			a(lenX-1)=-cofSec
			b(lenX-1)=1+2*cofSec+deltU(lenX-1)*deltT-cofFir*uOne(lenX-1)
		else
			a(lenX-1)=-cofSec-cofFir*uOne(lenX-1)
			b(lenX-1)=1+2*cofSec+deltU(lenX-1)*deltT+cofFir*uOne(lenX-1)
		end if

		do k=1,lenX-1
			r(k)=uOne(k)*deltU(k)*deltT+uZero(k)
		end do

		call ThomasAlgorithm(a,b,c,lenX-1,r)
		err=maxArr(r-uOne(1:lenX-1),lenX-1)	
		uOne(1:lenX-1)=r
		uOne(0)=bounLeft(time)                               ! Dirichlet boundary value
		uOne(lenX)=bounRight(time)		
	end do 

end subroutine SchemeSeven

!============================================================================================

subroutine SchemeEight(uOne,uZero,lenX,deltX,deltT,coeff,time)
! Crank-Nicolson scheme on page 347 (6.10.8)
! Newton's Method
! coeff is with the second order term
	use comNumFun
	real uOne(0:lenX),uZero(0:lenX)
	integer k
	real :: cofFir,cofSec,err         ! first order and second order's coeff in the diffition eqution 
	real a(2:lenX-1),b(1:lenX-1),c(1:lenX-2),r(1:lenX-1)
	
	err=1
	cofSec=coeff*deltT/deltX**2  
	cofFir=deltT/deltX
    uOne(0:lenX)=uZero(0:lenX)
    do while( abs(err)>1e-6 ) 
		b(1)=1+cofSec
		c(1)=-cofFir/4.*uOne(1)-cofSec/2.

		do k=2,lenX-2
			a(k)=cofFir/4.*uOne(k)-cofSec/2.
			b(k)=1+cofSec
			c(k)=-cofFir/4.*uOne(k)-cofSec/2.
		end do

		a(lenX-1)=cofFir/4.*uOne(lenX-1)-cofSec/2.
		b(lenX-1)=1+cofSec

		do k=1,lenX-1
			r(k)=-uOne(k)*cofFir/4.*(uOne(k+1)-uOne(k-1))+cofSec/2.*(uOne(k+1)-2*uOne(k)+uOne(k-1))+uZero(k)
		end do

		call ThomasAlgorithm(a,b,c,lenX-1,r)
		err=maxArr(r-uOne(1:lenX-1),lenX-1)	
		uOne(1:lenX-1)=r
		uOne(0)=bounLeft(time)                               ! Dirichlet boundary value
		uOne(lenX)=bounRight(time)		
	end do 

end subroutine SchemeEight