USMANTAURA
Programmer
program pisarchik
implicit none
real*8 t,tstep,y(5)
integer N,i,j,nsteps
N=5;nsteps=0.1;y(1)=0.6;y(2)=0.0;y(3)=0.6;y(4)=0.0;y(5)=0.0
open(6,file='p.dat')
write(6,*)0,Y(1),Y(2),Y(3),Y(4),Y(5)
do J = 1,nsteps
t = j*tstep
call rk4(t,y,tstep,N)
write(6,*)0,Y(1),Y(2),Y(3),Y(4),Y(5)
enddo
stop
end
subroutine rk4(t,y,tstep,N)
implicit none
real*8 deriv,h,t,tstep,y(5)
real*8 k1(5),k2(5),k3(5),k4(5),temp1(5),temp2(5),temp3(5)
real temp4(5),temp5(5)
integer i,N
h = tstep/2.0
do i = 1,N
k1(i) =tstep*deriv(t,y,i)
temp1(i) = y(i) + 0.5*k1(i)
enddo
do 20 i = 1,N
k2(i) =tstep*deriv(t+h,temp1,i)
temp2(i) = y(i) + 0.5*k2(i)
20 continue
do 30 i = 1,N
k2(i) =tstep*deriv(t+h,temp2,i)
temp3(i) = y(i) + k3(i)
30 continue
do 40 i = 1,N
k4(i) =tstep*deriv(t+tstep,temp3,i)
y(i) = y(i) + k1(i) + 2.0*(k2(i) + k3(i) + k4(i))/6.0
40 continue
return
end
function(t,temp,i)
implicit none
real*8 deriv,t,temp(5)
integer i
if(i.eq.1)deriv = temp(2)
if(i.eq.2)deriv = -0.4*temp(2)+(0.25+temp(5)) - 0.5*temp(1)**3
& - delta*temp(1)*temp(3)**2
if(i.eq.3)deriv = temp(4)
if(i.eq.4)deriv = -0.4*temp(4)+0.25*temp(3) - 0.5*temp(3)**3
& - delta*(temp(1)**2)*temp(3)
if(i.eq.5)deriv = 2*3.14*sqrt((0.25*0.8)**2-temp(5)**2)
return
end
System to be solve as an example is
Want δ to vary from 0.3 to 0.5 and
1) plot x1 versus δ
2) x3 versus δ
implicit none
real*8 t,tstep,y(5)
integer N,i,j,nsteps
N=5;nsteps=0.1;y(1)=0.6;y(2)=0.0;y(3)=0.6;y(4)=0.0;y(5)=0.0
open(6,file='p.dat')
write(6,*)0,Y(1),Y(2),Y(3),Y(4),Y(5)
do J = 1,nsteps
t = j*tstep
call rk4(t,y,tstep,N)
write(6,*)0,Y(1),Y(2),Y(3),Y(4),Y(5)
enddo
stop
end
subroutine rk4(t,y,tstep,N)
implicit none
real*8 deriv,h,t,tstep,y(5)
real*8 k1(5),k2(5),k3(5),k4(5),temp1(5),temp2(5),temp3(5)
real temp4(5),temp5(5)
integer i,N
h = tstep/2.0
do i = 1,N
k1(i) =tstep*deriv(t,y,i)
temp1(i) = y(i) + 0.5*k1(i)
enddo
do 20 i = 1,N
k2(i) =tstep*deriv(t+h,temp1,i)
temp2(i) = y(i) + 0.5*k2(i)
20 continue
do 30 i = 1,N
k2(i) =tstep*deriv(t+h,temp2,i)
temp3(i) = y(i) + k3(i)
30 continue
do 40 i = 1,N
k4(i) =tstep*deriv(t+tstep,temp3,i)
y(i) = y(i) + k1(i) + 2.0*(k2(i) + k3(i) + k4(i))/6.0
40 continue
return
end
function(t,temp,i)
implicit none
real*8 deriv,t,temp(5)
integer i
if(i.eq.1)deriv = temp(2)
if(i.eq.2)deriv = -0.4*temp(2)+(0.25+temp(5)) - 0.5*temp(1)**3
& - delta*temp(1)*temp(3)**2
if(i.eq.3)deriv = temp(4)
if(i.eq.4)deriv = -0.4*temp(4)+0.25*temp(3) - 0.5*temp(3)**3
& - delta*(temp(1)**2)*temp(3)
if(i.eq.5)deriv = 2*3.14*sqrt((0.25*0.8)**2-temp(5)**2)
return
end
System to be solve as an example is
Want δ to vary from 0.3 to 0.5 and
1) plot x1 versus δ
2) x3 versus δ