Απο το θεμα 8 που εκανε ο Αναγνωστοπουλος πρεπει να ειναι τα ιδια.
Για την δυναμικη ενεργεια U = - ολοκληρωμα (F*dx) ==> U/m= - ολοκληρωμα (a (x)) , οπου a(x) ειναι η επιταχυνση. Οποτε κανεις στο χαρτι το ολοκληρωμα κ το δινεις να το βαρεσει για t......dt το προγραμμα.
Ολικη ενεργεια = Δυναμικη + Κινητικη ==> F/m = T/m + U/m = 1/2* U^2 + U όπου η κινητική ενεργεια ειναι 1/2*U^2 δλδ το μισο του τετραγωνου της δυναμικης..
Σου δειχνω και το προγραμμα του θεματος 8 για να δεις πως βαζει τις ενεργειες. Με παρομοιο τροπο μπαινουν και στο θεμα 7
Κώδικας: Επιλογή όλων
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Program to integrate equations of motion for accelerations
C which are functions of x with the method of Euler, Euler-Cromer
C and Euler-Verlet.
C The user sets initial conditions and the subroutines return
C X(t) and V(t)=dX(t)/dt in arrays T(1..STEPS),X(1..STEPS),V(1..STEPS)
C The user provides number of integration STEPS and the final
C time TFI.Initial time is assumed to be t_0=0 and the integration
C step h = TFI/(STEPS-1)
C The user programs a real function accel(x) which gives the
C acceleration dV(t)/dt as function of X.
C NOTE: T(1) = 0 T(STEPS) = TFI and there are STEPS-1 aditional
C steps after the initial point
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
program diff_eq_euler
implicit none ! We force ourselved to declare all variables
integer P ! The size of the arrays, should be larger
parameter(P=110000) ! than number of steps
real*8 T(P),X(P),V(P) ! time t,x(t),v(t)=dx/dt
real*8 Xin,Vin,Tfi ! initial conditions
real*8 Energy
real*8 KineticE,PotentialE,accel !Allages!!
integer steps,i
!t_0 = 0.0
C The user provides initial conditions X_0,V_0 final time t_f and
C number of steps:
print *,'Enter X_0,V_0,t_f,number of steps (t_0=0):'
read(5,*)Xin,Vin,Tfi,steps
C This check is necessary to avoid memory violations:
if(steps .ge. p )then
print *,'steps must be strictly less than p. steps,p= ',steps,p
stop
endif
C Xin= X(1), Vin=V(1), T(1)=0 and the routine gives eveolution in
C T(2..STEPS), X(2..STEPS), V(2..STEPS) which we print in a file
call euler(Xin,Vin,Tfi,steps,T,X,V)
open(unit=20,file="euler.dat") !filename euler.dat given here
do i=1,steps
C Each line in data file has time, position, velocity:
Energy = 0.5D0*V(i)*V(i)
$ + 0.5D0*39.47841760D0*X(i)*X(i)
$ + 0.25D0*20.0D0*X(i)*X(i)*X(i)*X(i)
$ + 10.0D0*X(i)
write(20,*) T(i),X(i),V(i),Energy
enddo
close(20) !we close the unit to be reused below
C We repeat everything for each method
call euler_cromer(Xin,Vin,Tfi,steps,T,X,V)
open(unit=20,file="euler_cromer.dat")
do i=1,steps
Energy = 0.5D0*V(i)*V(i)
$ + 0.5D0*39.47841760D0*X(i)*X(i)
$ + 0.25D0*20.0D0*X(i)*X(i)*X(i)*X(i)
$ + 10.0D0*X(i)
write(20,*) T(i),X(i),V(i),Energy
enddo
close(20)
call euler_verlet(Xin,Vin,Tfi,steps,T,X,V)
open(unit=20,file="euler_verlet.dat")
open(unit=21,file="energy.dat")
do i=1,steps
KineticE = 0.5D0*V(i)*V(i)
PotentialE = 0.5D0*39.47841760D0*X(i)*X(i)
$ + 0.25D0*20.0D0*X(i)*X(i)*X(i)*X(i)
$ + 10.0D0*X(i)
Energy = KineticE + PotentialE
write(20,*) T(i),X(i),V(i),Energy
write(21,*) T(i),X(i),accel(X(i)),KineticE,PotentialE,Energy
enddo
close(20)
close(21)
end
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Function which returns the value of acceleration at
C position x used in the integration subroutines
C euler, euler_cromer and euler_verlet
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
real*8 function accel(x)
implicit none
real*8 x
accel = -39.47841760D0*x-20.0D0*x*x*x-10.0D0
end
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Driver routine for integrating equations of motion
C using the Euler method
C Input:
C Xin=X(1), Vin=V(1) -- initial condition at t=0,
C Tfi the final time and steps the number of steps of integration
C (the initial point is counted as the first one)
C Output:
C The arrays T(1..steps), X(1..steps), V(1..steps) which
C gives x(t_i)=X(i), dx/dt(t_i)=V(i), t_i=T(i) i=1..steps
C where for i=1 we have the initial condition.
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine euler(Xin,Vin,Tfi,steps,T,X,V)
implicit none
integer P
parameter(P=110000)
real*8 T(P),X(P),V(P) !time t,x(t),v(t)=dx/dt
real*8 Xin,Vin,Tfi
integer steps,i
real*8 h,accel !**Note: we have to declare the function accel**
C Initial conditions set here:
T(1) = 0.0D0
X(1) = Xin
V(1) = Vin
C h is the time step Dt
h = Tfi/(steps-1)
do i = 2,steps
T(i) = T(i-1)+h ! time advances by Dt
X(i) = X(i-1)+V(i-1)*h ! advancement and storage of position
V(i) = V(i-1)+accel(X(i-1))*h !... and velocity. Here we call accel(x)
enddo
end
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Driver routine for integrating equations of motion
C using the Euler-Cromer method
C Input:
C Xin=X(1), Vin=V(1) -- initial condition at t=0,
C Tfi the final time and steps the number of steps of integration
C (the initial point is counted as the first one)
C Output:
C The arrays T(1..steps), X(1..steps), V(1..steps) which
C gives x(t_i)=X(i), dx/dt(t_i)=V(i), t_i=T(i) i=1..steps
C where for i=1 we have the initial condition.
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine euler_cromer(Xin,Vin,Tfi,steps,T,X,V)
implicit none
integer P
parameter(P=110000)
real*8 T(P),X(P),V(P) !time t,x(t),v(t)=dx/dt
real*8 Xin,Vin,Tfi
integer steps,i
real*8 h,accel
T(1) = 0.0
X(1) = Xin
V(1) = Vin
h = Tfi/(steps-1)
do i = 2,steps
T(i) = T(i-1)+h
V(i) = V(i-1)+accel(X(i-1))*h
X(i) = X(i-1)+V(i)*h !here is the difference compared to Euler
enddo
end
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Driver routine for integrating equations of motion
C using the Euler - Verlet method
C Input:
C Xin=X(1), Vin=V(1) -- initial condition at t=0,
C Tfi the final time and steps the number of steps of integration
C (the initial point is counted as the first one)
C Output:
C The arrays T(1..steps), X(1..steps), V(1..steps) which
C gives x(t_i)=X(i), dx/dt(t_i)=V(i), t_i=T(i) i=1..steps
C where for i=1 we have the initial condition.
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine euler_verlet(Xin,Vin,Tfi,steps,T,X,V)
implicit none
integer P
parameter(P=110000)
real*8 T(P),X(P),V(P) !time t,x(t),v(t)=dx/dt
real*8 Xin,Vin,Tfi
integer steps,i
real*8 h,g_over_l
parameter(g_over_l=10.0D0)
real*8 h2,X0,o2h
real*8 accel
C Initial conditions set here:
T(1) = 0.0D0
X(1) = Xin
V(1) = Vin
h = Tfi/(steps-1) ! time step
h2 = h*h ! time step squared
o2h = 0.5D0/h ! h/2
C We have to initialize one more step: X0 corresponds to 'X(0)'
X0 = X(1) - V(1) * h + accel(X(1)) *h2/2.0D0
T(2) = h
X(2) = 2.0D0*X(1) - X0 + accel(X(1)) *h2
C Now i starts from 3:
do i = 3,steps
T(i) = T(i-1)+h
X(i) = 2.0D0*X(i-1) - X(i-2) + accel(X(i-1))*h2
V(i-1) = o2h * (X(i)-X(i-2))
enddo
C Notice that we have one more step for the velocity:
V(steps)= (X(steps)-X(steps-1))/h
end
Facebook
Reddit
Discord Invitation
