Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/difequ.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

difequ.f90
View Options

!>
!> @file difequ.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine difequ ( mode, xx, v )
!*************************************************************************
!
!! DIFEQU returns information about a differential equation.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, integer MODE, an integer indicating the task to be performed.
! 1, initialization
! 2, evaluate de at xx
! 3, specify the next side condition
! 4, analyze the approximation
!
! Input, real ( kind = 8 ) XX, a point at which information is wanted
!
! Output, real ( kind = 8 ) V, depends on the MODE. See comments below
!
implicit none
integer, parameter :: npiece = 100
integer, parameter :: ncoef = 2000
real ( kind = 8 ) break
real ( kind = 8 ) coef
real ( kind = 8 ) eps
real ( kind = 8 ) ep1
real ( kind = 8 ) ep2
real ( kind = 8 ) error
real ( kind = 8 ) factor
integer i
integer iside
integer itermx
integer k
integer kpm
integer l
integer m
integer mode
real ( kind = 8 ) rho
real ( kind = 8 ) s2ovep
real ( kind = 8 ) solutn
real ( kind = 8 ) un
real ( kind = 8 ) v(20)
real ( kind = 8 ) value
real ( kind = 8 ) x
real ( kind = 8 ) xside
real ( kind = 8 ) xx
common /approx/ break(npiece),coef(ncoef),l,kpm
common /side/ m,iside,xside(10)
common /other/ itermx,k,rho(19)
!
! This sample of DIFEQU is for the example in chapter xv. It is a
! nonlinear second order two point boundary value problem.
!
go to (10,50,60,110), mode
!
! Initialize everything, Set the order M of the differential equation,
! the nondecreasing sequence xside(i),i=1,...,m, of points at which side
! conditions are given and anything else necessary.
!
10 continue
m = 2
xside(1) = 0.0D+00
xside(2) = 1.0D+00
!
! Print out heading.
!
write ( *, * ) ' '
write ( *, * ) ' Carrier''s nonlinear perturb. problem'
write ( *, * ) ' '
eps = 0.005D+00
write(*,*)'EPS = ',eps
!
! Set constants used in formula for solution below.
!
factor = ( sqrt ( 2.0D+00 ) + sqrt ( 3.0D+00 ) )**2
s2ovep = sqrt ( 2.0D+00 / eps )
!
! Initial guess for Newton iteration. un(x)=x*x-1.
!
l = 1
break(1) = 0.0D+00
do i = 1, kpm
coef(i) = 0.0D+00
end do
coef(1) = -1.0D+00
coef(3) = 2.0D+00
itermx = 10
return
!
! Provide value of left side coefficients and right side at xx .
! specifically, at xx the dif.equ. reads:
!
! v(m+1)d**m+v(m)d**(m-1) + ... + v(1)d**0 = v(m+2)
!
! in terms of the quantities v(i),i=1,...,m+2, to be computed here.
!
50 continue
v(3) = eps
v(2) = 0.0D+00
call ppvalu(break,coef,l,kpm,xx,0,un)
v(1) = 2.0D+00 * un
v(4) = un**2 + 1.0D+00
return
!
! provide the M side conditions. these conditions are of the form
! v(m+1)d**m+v(m)d**(m-1) + ... + v(1)d**0 = v(m+2)
! in terms of the quantities v(i),i=1,...,m+2, to be specified here.
! note that v(m+1)=0 for customary side conditions.
!
60 continue
v(m+1) = 0.0D+00
if ( iside == 1 ) then
v(2) = 1.0D+00
v(1) = 0.0D+00
v(4) = 0.0D+00
iside = iside+1
else if ( iside == 2 ) then
v(2) = 0.0D+00
v(1) = 1.0D+00
v(4) = 0.0D+00
iside = iside + 1
end if
return
!
! Calculate the error near the boundary layer at 1.
!
110 continue
write(*,*)' '
write(*,*)' X, G(X) and G(X)-F(X) at selected points:'
write(*,*)' '
x = 0.75D+00
do i = 1, 9
ep1 = exp ( s2ovep * ( 1.0D+00 - x ) ) * factor
ep2 = exp ( s2ovep * ( 1.0D+00 + x ) ) * factor
solutn = 12.0D+00 / ( 1.0D+00 + ep1 )**2 * ep1 &
+12.0D+00 / ( 1.0D+00 + ep2 )**2 * ep2 - 1.0D+00
call ppvalu(break,coef,l,kpm,x,0,value)
error = solutn-value
write ( *, '(1x,3g14.6)' ) x, solutn, error
x = x+0.03125
end do
return
end
c4science · Help