need help

[rasha80]

Am using Fortran to solve numerical equations using Euler and Runge-Kutta methods and need an expert opinion in the programme I did, can anyone help me plz??

the euler programme

PROGRAM membrane
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,con,edon
DOUBLE PRECISION deex,xstep,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2),con(10,2),edon(10,2)
COMMON/a/c(20001,2),pot(20001),conf(20001,2)
COMMON/beth/farad,gascon,deex,umem,nep,xstep,temp,node

OPEN(unit=5,file='memb.dat',status='unknown')
OPEN(unit=6,file='memb.res',status='unknown')
farad=96485.3383
gascon=8.3145
! convrg=1.0D-5
! PRINT*,' auto 1 other 0'
! read*,iauto
iauto = 1
IF(iauto.EQ.1)THEN
READ(5,*)nep
print*, nep
write(6,*) nep
READ(5,*)temp,umem
print*,temp,umem
write(6,*)temp,umem
DO i = 1,nep
DO j = 1,2
READ(5,*)vpm(i,j),diff(i,j),hidc(i,j),edon(i,j),con(i,j)
PRINT*,i,j,vpm(i,j),diff(i,j),hidc(i,j),edon(i,j),con(i,j)
WRITE(6,*)i,j,vpm(i,j),diff(i,j),hidc(i,j)
WRITE(6,*)i,j,edon(i,j),con(i,j)
END DO
vpm(i,1)=DABS(vpm(i,1))
vpm(i,2)=-DABS(vpm(i,2))
END DO
ELSE
Print*,' Where is the data??'
stop
END IF
node=2000
! thick = membrane thickness in micron
deex = 10.0
deex = deex*1.0D-6
xstep = deex/(DFLOAT(node))
print*,' '
print*,'**************************************************'
print*,' '
print*,' ALL UNITS ARE S.I. '
print*,' '
print*,'**************************************************'
print*,' '
! set initial conditions
pot(1) = 0.0
CALL euler
CALL output
STOP
END

SUBROUTINE EULER
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION conf,vpm,diff,hidc,temp,umem,con,edon
DOUBLE PRECISION deex,xstep,farad,gascon
DOUBLE PRECISION change,dipsidxn,dipsidxnol
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2),con(10,2),edon(10,2)
COMMON/a/c(20001,2),pot(20001),conf(20001,2)
COMMON/beth/farad,gascon,deex,umem,nep,xstep,temp,node
dipsidxn = 0.0
DO i=2,node+1
change = 1.0
! set initial conditions
conf(1,1)=con(1,1)
conf(1,2)=con(1,2)
while(change.GT.0.001)DO
c(1,1)=con(1,1)*exp(-vpm(1,1)*farad*edon(1,1)/gascon/temp)
c(1,2)=con(1,2)*exp(-vpm(1,2)*farad*edon(1,2)/gascon/temp)
const1 = (umem/diff(1,1))*(hidc(1,1)*c(i-1,1)-conf(i-1,1))
const3 = (umem/diff(1,2))*(hidc(1,2)*c(i-1,2)-conf(i-1,2))
const2 = (vpm(1,1)*c(i-1,1)*farad/gascon/temp)*dipsidxn
const4 = (vpm(1,2)*c(i-1,2)*farad/gascon/temp)*dipsidxn
c(i,1) = c(i-1,1) + xstep*(const1-const2)
c(i,2) = c(i-1,2) + xstep*(const3-const4)
sum1 = 0.0
sum2 = 0.0
DO k=1,2
DO j=1,nep
const5 = vpm(j,k)*umem/diff(j,k)
const6 = hidc(j,k)*c(i-1,j)-conf(j,k)
sum1 = sum1 + const5*const6
sum2 = sum2+vpm(j,k)*vpm(j,k)*c(i-1,j)
END DO
END DO
sum2 = sum2*farad/gascon/temp
pot(i) = pot(i-1) + xstep*sum1/sum2
dipsidxn = sum1/sum2
change = Dabs((dipsidxn - dipsidxnol)/dipsidxn)
dipsidxnol = dipsidxn
cont5=(vpm(1,1)*c(i-1,1)*farad/gascon/temp)*(sum1/sum2)
cont6=hidc(1,1)*c(i-1,1)
cont7=diff(1,1)/umem
conf(i,1)=(cont6)-((((c(i,1)-c(i-1,1))/xstep)+(cont5))*cont7)
cont8=(vpm(1,2)*c(i-1,2)*farad/gascon/temp)*(sum1/sum2)
cont9=hidc(1,2)*c(i-1,2)
cont10=diff(1,2)/umem
conf(i,2)=(cont9)-((((c(i,2)-c(i-1,2))/xstep)+(cont8))*cont10)
IF(change.lt.0.001)then
print*, i,change
end if
END DO
END DO
RETURN
END

SUBROUTINE output
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION conf,vpm,diff,hidc,temp,umem,con,edon
DOUBLE PRECISION deex,xstep,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2),con(10,2),edon(10,2)
COMMON/a/c(20001,2),pot(20001),conf(20001,2)
COMMON/beth/farad,gascon,deex,umem,nep,xstep,temp,node
i = 1
write(6,10) i,c(1,1),c(1,2),pot(1),conf(1,1),conf(1,2)
DO i = 2,node+1
iz = i/40
ii = 40*iz
IF(i.eq.ii)THEN
print*, i,c(i,1),c(i,2),pot(i),conf(i,1),conf(i,2)
write(6,10) i,c(i,1),c(i,2),pot(i),conf(i,1),conf(i,2)
END IF
END DO
10 FORMAT(i6,5(5x,e13.5))
RETURN
END

 

[xwb]

What kind of opinion: method used, code layout, coding technique?

 

[mikrom]

Hi rasha80,

On the programming style:
IMHO: the COMMON blocks are obsolete. It would be ok when you modified an old program, but when I should write a new one I could never use COMMON blocks for passing arguments to procedures/functions because IMHO it's not good style and maintenace of such programms is terrible.
Rather look at the newer Fortran dialects and pass the arguments via (..) and intent(in)/(intent(out).

On the numerics:
I have no tried your program, so excuse me if I'm wrong. But it seems for me, that you are using constant step. If you want integrate the ordinary differential equation properly, you should use step control. Because on some subintervals you can integrate it with larger step, but sometimes there are subintervalls (especially by stiff equations) where you need to decrease your step to get reasonable results.

 

[rasha80]

thanks for your posts, i enclosed the two programmes with more details about the work. i really need help and an expert opinion in this work. is there a way to get in contact with you in order to explain it in more details???? thanks alot for your help.






8.1.1 Euler’s numerical method
At first equation (8.3) was written according to Euler’s method as follows
(8.1)
From equation (8.17) the concentration inside the membrane and the permeate concentration can be calculated. Then equation (8.9) was rearranged and used to determine the initial solute concentration inside the membrane (ci) by using the solute feed concentration. Equation (8.9) was rearranged and written as
(8.2)
where Ci,f is the solute concentration in the feed solution. Then equation (8.6) was used to calculate the potential gradient (d?/dx), where it was substituted into equation (8.17) to calculate a new value for the solute concentration inside the membrane. The new value of the solute concentration inside the membrane was used to calculate the solute permeate concentration, by substituting it into equation (8.17). The initial permeate concentration was assumed to be equal to the feed concentration to be able to calculate the initial value of the concentration inside the membrane. The program used (program-1) is given in appendix-1.
Euler’s method was also used to calculate the permeate concentration and the concentration inside the membrane by using equation (8.16) instead of equation (8.3). Where equation (8.16) was written according to Euler’s method as follows
(8.3)
The initial permeate concentration was assumed to be equal to the feed concentration. Then equation (8.18) was used to determine the initial solute concentration inside the membrane. Then equation (8.6) was used to calculate the potential gradient (d?/dx), where it was substituted into equation (8.19) to calculate a new value for the solute concentration inside the membrane. The new value of the solute concentration inside the membrane was used to calculate the solute permeate concentration, by substituting it into equation (8.19). The program used (program-2) is given in appendix-2.
8.1.2 Runge-Kutta numerical method
The Runge-Kutta method was used to calculate the concentration change inside the membrane. At first the initial concentration inside the membrane was assumed to be equal to the permeate concentration. Then the potential gradient was calculated by using equation (8.6). Then the potential gradient was substituted into equation (8.3) where it was integrated and a new value of concentration inside the membrane resulted. The program used (program-3) is given in appendix-3.
(8.4)
(8.5)


The Fortran programmes are enclosed in below where runge-kutta is runge-kutta1 and runge-kutta2, and appendix 1 is for Euler method. The difference between runge-kutta 1 and 2 that conf is recalculated rather than being considered constant (but this programme is not working were it start but stops and do not know where the error is, can someone help me how to find the error please).

My question is there a way to write the equations in Euler programme as in runge-kutta?

c(1,1)=con(1,1)*exp(-vpm(1,1)*farad*edon(1,1)/gascon/temp)
c(1,2)=con(1,2)*exp(-vpm(1,2)*farad*edon(1,2)/gascon/temp)
const1 = (umem/diff(1,1))*(hidc(1,1)*c(i-1,1)-conf(i-1,1))
const3 = (umem/diff(1,2))*(hidc(1,2)*c(i-1,2)-conf(i-1,2))
const2 = (vpm(1,1)*c(i-1,1)*farad/gascon/temp)*dipsidxn
const4 = (vpm(1,2)*c(i-1,2)*farad/gascon/temp)*dipsidxn
c(i,1) = c(i-1,1) + xstep*(const1-const2)
c(i,2) = c(i-1,2) + xstep*(const3-const4)
sum1 = 0.0
sum2 = 0.0
DO k=1,2
DO j=1,nep
const5 = vpm(j,k)*umem/diff(j,k)
const6 = hidc(j,k)*c(i-1,k)-conf(i-1,k)
sum1 = sum1 + const5*const6
sum2 = sum2+vpm(j,k)*vpm(j,k)*c(i-1,k)
END DO
END DO
sum2 = sum2*farad/gascon/temp
pot(i) = pot(i-1) + xstep*sum1/sum2
dipsidxn = sum1/sum2
dipsidxnol = dipsidxn
change = Dabs((dipsidxn - dipsidxnol)/dipsidxn)
cont5=(vpm(1,1)*c(i-1,1)*farad/gascon/temp)*(sum1/sum2)
cont6=hidc(1,1)*c(i-1,1)
cont7=diff(1,1)/umem
conf(i,1)=(cont6)-((((c(i,1)-c(i-1,1))/xstep)+(cont5))*cont7)
cont8=(vpm(1,2)*c(i-1,2)*farad/gascon/temp)*(sum1/sum2)
cont9=hidc(1,2)*c(i-1,2)
cont10=diff(1,2)/umem
conf(i,2)=(cont9)-((((c(i,2)-c(i-1,2))/xstep)+(cont8))*cont10)


Where the red colour lines is equation 8.6. Where it must be solved in order to solve equation 8.3.




Euler method:


PROGRAM membrane
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,con,edon
DOUBLE PRECISION deex,xstep,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2),con(10,2),edon(10,2)
COMMON/a/c(20001,2),pot(20001),conf(20001,2)
COMMON/beth/farad,gascon,deex,umem,nep,xstep,temp,node

OPEN(unit=5,file='memb.dat',status='unknown')
OPEN(unit=6,file='memb.res',status='unknown')
farad=96485.3383
gascon=8.3145
! convrg=1.0D-5
! PRINT*,' auto 1 other 0'
! read*,iauto
iauto = 1
IF(iauto.EQ.1)THEN
READ(5,*)nep
print*, nep
write(6,*) nep
READ(5,*)temp,umem
print*,temp,umem
write(6,*)temp,umem
DO i = 1,nep
DO j = 1,2
READ(5,*)vpm(i,j),diff(i,j),hidc(i,j),edon(i,j),con(i,j)
PRINT*,i,j,vpm(i,j),diff(i,j),hidc(i,j),edon(i,j),con(i,j)
WRITE(6,*)i,j,vpm(i,j),diff(i,j),hidc(i,j)
WRITE(6,*)i,j,edon(i,j),con(i,j)
END DO
vpm(i,1)=DABS(vpm(i,1))
vpm(i,2)=-DABS(vpm(i,2))
END DO
ELSE
Print*,' Where is the data??'
stop
END IF
node=200
! thick = membrane thickness in micron
deex = 10.0
deex = deex*1.0D-6
xstep = deex/(DFLOAT(node))
print*,' '
print*,'**************************************************'
print*,' '
print*,' ALL UNITS ARE S.I. '
print*,' '
print*,'**************************************************'
print*,' '
! set initial conditions
pot(1) = 0.0
CALL euler
CALL output
STOP
END

SUBROUTINE EULER
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION conf,vpm,diff,hidc,temp,umem,con,edon
DOUBLE PRECISION deex,xstep,farad,gascon
DOUBLE PRECISION change,dipsidxn,dipsidxnol
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2),con(10,2),edon(10,2)
COMMON/a/c(20001,2),pot(20001),conf(20001,2)
COMMON/beth/farad,gascon,deex,umem,nep,xstep,temp,node
dipsidxn = 0.0
DO i=2,node+1
change = 1.0
! set initial conditions
conf(1,1)=con(1,1)
conf(1,2)=con(1,2)
while(change.GT.0.001)DO
c(1,1)=con(1,1)*exp(-vpm(1,1)*farad*edon(1,1)/gascon/temp)
c(1,2)=con(1,2)*exp(-vpm(1,2)*farad*edon(1,2)/gascon/temp)
const1 = (umem/diff(1,1))*(hidc(1,1)*c(i-1,1)-conf(i-1,1))
const3 = (umem/diff(1,2))*(hidc(1,2)*c(i-1,2)-conf(i-1,2))
const2 = (vpm(1,1)*c(i-1,1)*farad/gascon/temp)*dipsidxn
const4 = (vpm(1,2)*c(i-1,2)*farad/gascon/temp)*dipsidxn
c(i,1) = c(i-1,1) + xstep*(const1-const2)
c(i,2) = c(i-1,2) + xstep*(const3-const4)
sum1 = 0.0
sum2 = 0.0
DO k=1,2
DO j=1,nep
const5 = vpm(j,k)*umem/diff(j,k)
const6 = hidc(j,k)*c(i-1,k)-conf(i-1,k)
sum1 = sum1 + const5*const6
sum2 = sum2+vpm(j,k)*vpm(j,k)*c(i-1,k)
END DO
END DO
sum2 = sum2*farad/gascon/temp
pot(i) = pot(i-1) + xstep*sum1/sum2
dipsidxn = sum1/sum2
dipsidxnol = dipsidxn
change = Dabs((dipsidxn - dipsidxnol)/dipsidxn)
cont5=(vpm(1,1)*c(i-1,1)*farad/gascon/temp)*(sum1/sum2)
cont6=hidc(1,1)*c(i-1,1)
cont7=diff(1,1)/umem
conf(i,1)=(cont6)-((((c(i,1)-c(i-1,1))/xstep)+(cont5))*cont7)
cont8=(vpm(1,2)*c(i-1,2)*farad/gascon/temp)*(sum1/sum2)
cont9=hidc(1,2)*c(i-1,2)
cont10=diff(1,2)/umem
conf(i,2)=(cont9)-((((c(i,2)-c(i-1,2))/xstep)+(cont8))*cont10)
IF(change.lt.0.001)then
print*, i,change
end if
END DO
END DO
RETURN
END

SUBROUTINE output
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION conf,vpm,diff,hidc,temp,umem,con,edon
DOUBLE PRECISION deex,xstep,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2),con(10,2),edon(10,2)
COMMON/a/c(20001,2),pot(20001),conf(20001,2)
COMMON/beth/farad,gascon,deex,umem,nep,xstep,temp,node
i = 1
write(6,10) i,c(1,1),c(1,2),pot(1),conf(1,1),conf(1,2)
DO i = 2,node+1
iz = i/40
ii = 40*iz
IF(i.eq.ii)THEN
print*, i,c(i,1),c(i,2),pot(i),conf(i,1),conf(i,2)
write(6,10) i,c(i,1),c(i,2),pot(i),conf(i,1),conf(i,2)
END IF
END DO
10 FORMAT(i6,5(5x,e13.5))
RETURN
END


Runge-Kutta:

PROGRAM membrane runge kutta
C driver for routine membrane
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2)
COMMON/a/conf(20001,2)
COMMON/beth/farad,gascon,umem,nep,temp,node
INTEGER NSTEP,NVAR
PARAMETER(NVAR=4,NSTEP=100)
INTEGER i,j
REAL xx(200),x1,x2,y(2,200),vstart(NVAR)
COMMON /path/ xx,y
EXTERNAL derivs

OPEN(unit=5,file='memb.dat',status='unknown')
OPEN(unit=6,file='memb.res',status='unknown')

farad=96485.3383
gascon=8.3145
! PRINT*,' auto 1 other 0'
! read*,iauto
iauto = 1
IF(iauto.EQ.1)THEN
READ(5,*)nep
print*, nep
write(6,*) nep
READ(5,*)temp,umem
print*,temp,umem
write(6,*)temp,umem
DO i = 1,nep
DO j = 1,2
READ(5,*)vpm(i,j),diff(i,j),hidc(i,j)
PRINT*,i,j,vpm(i,j),diff(i,j),hidc(i,j)
WRITE(6,*)i,j,vpm(i,j),diff(i,j),hidc(i,j)
END DO
vpm(i,1)=DABS(vpm(i,1))
vpm(i,2)=-DABS(vpm(i,2))
END DO
ELSE
Print*,' Where is the data??'
stop
END IF
node=100
y(1,1)=3.0D-3
y(2,1)=3.0D-3
conf(1,1)=1.0D-3
conf(1,2)=1.0D-3
! thick = membrane thickness in micron
x1=0.0D-6
vstart(1)=y(1,1)
vstart(2)=y(2,1)
x2=20.0D-6
! set intial value of y
call membrane(vstart,NVAR,x1,x2,NSTEP,derivs)
call output
do 11 i=1,100
print*,i,xx(i),y(1,i),y(2,i),conf(i,1),conf(i,2)
11 continue
END

SUBROUTINE derivs(x,y,dydx)
REAL x,y(*),dydx(*)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2)
COMMON/a/conf(20001,2)
COMMON/beth/farad,gascon,umem,nep,temp,node
DO i=1,node+1
change = 1.0
while(change.GT.0.001)DO
const1=(vpm(1,1)*umem/diff(1,1))*(hidc(1,1)*y(1)-conf(i,1))
const2=(vpm(1,2)*umem/diff(1,2))*(hidc(1,2)*y(2)-conf(i,2))
const3=vpm(1,1)*vpm(1,1)*y(1)
const4=vpm(1,2)*vpm(1,2)*y(2)
cons=(const1+const2)/((farad/gascon/temp)*(const3+const4))
const5=(umem/diff(1,1))*(hidc(1,1)*y(1)-conf(i,1))
const6=(vpm(1,1)*y(1)*farad/gascon/temp)*cons
const7=(umem/diff(1,2))*(hidc(1,2)*y(2)-conf(i,2))
const8=(vpm(1,2)*y(2)*farad/gascon/temp)*cons
dydx(1)=(const5-const6)
dydx(2)=(const7-const8)
cont1=(vpm(1,1)*y(1)*farad/gascon/temp)*cons
cont2=hidc(1,1)*y(1)
cont3=diff(1,1)/umem
conf(i+1,1)=(cont2)-((dydx(1)+cont1)*cont3)
cont4=(vpm(1,2)*y(2)*farad/gascon/temp)*cons
cont5=hidc(1,2)*y(2)
cont6=diff(1,2)/umem
conf(i+1,2)=(cont5)-((dydx(2)+cont4)*cont16)
IF(change.lt.0.001)then
print*, i,change
end if
END DO
END DO
return
END

SUBROUTINE membrane(vstart,nvar,x1,x2,nstep,derivs)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2)
COMMON/a/conf(20001,2)
COMMON/beth/farad,gascon,umem,nep,temp,node
INTEGER nstep,nvar,NMAX,NSTPMX
PARAMETER (NMAX=4,NSTPMX=200)
REAL x1,x2,vstart(nvar),xx(NSTPMX),y(NMAX,NSTPMX)
EXTERNAL derivs
COMMON /path/ xx,y
CU USES rk4
INTEGER i,k
REAL h,x,dv(NMAX),v(NMAX)
do 11 i=1,nvar
v(i)=vstart(i)
y(i,1)=v(i)
11 continue
xx(1)=x1
x=x1
h=(x2-x1)/nstep
do 13 k=1,nstep
call derivs(x,v,dv)
call rk4(v,dv,nvar,x,h,v,derivs)
if(x+h.eq.x)pause 'stepsize not significant in rkdumb'
x=x+h
xx(k+1)=x
do 12 i=1,nvar
y(i,k+1)=v(i)
12 continue
print*,xx(k+1),y(1,k+1),y(2,k+1)
print*,'h',h
print*,'nstep',nstep
13 continue
return
END

SUBROUTINE rk4(y,dydx,n,x,h,yout,derivs)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2)
COMMON/a/conf(20001,2)
COMMON/beth/farad,gascon,umem,nep,temp,node
INTEGER n,NMAX
REAL h,x,dydx,y,yout
EXTERNAL derivs
PARAMETER (NMAX=4)
INTEGER i
REAL h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX)
hh=h*0.5
h6=h/6.
xh=x+hh
do 11 i=1,n
yt(i)=y(i)+hh*dydx(i)
11 continue
call derivs(xh,yt,dyt)
do 12 i=1,n
yt(i)=y(i)+hh*dyt(i)
12 continue
call derivs(xh,yt,dym)
do 13 i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
13 continue
call derivs(x+h,yt,dyt)
do 14 i=1,n
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
14 continue
return
END

SUBROUTINE output
IMPLICIT REAL (A-H,O-Z)
DOUBLE PRECISION vpm,diff,hidc,temp,umem,farad,gascon
COMMON/aleph/vpm(10,2),diff(10,2),hidc(10,2)
COMMON/a/conf(20001,2)
COMMON/beth/farad,gascon,umem,nep,temp,node
INTEGER NSTEP,NVAR
PARAMETER(NVAR=4,NSTEP=20)
INTEGER i,j
REAL xx(200),x1,x2,y(2,200),vstart(NVAR)
COMMON /path/ xx,y
DO i = 1,node+1
iz = i/1
ii = 1*iz
IF(i.eq.ii)THEN
print*, i,xx(i),y(1,i),y(2,i),conf(i,1),conf(i,2)
write(6,10)i,xx(i),y(1,i),y(2,i),conf(i,1),conf(i,1)
END IF
END DO
10 FORMAT(i6,5(5x,e13.7))
RETURN
END