МЕХАНИЧЕСКИЙ ЯЗЫК ПРОГРАММИРОВАНИЯ
дифференциальные уравнения только для обучения, сделанные
дифференциальные уравнения только для обучения, сделанные
C:
module precision
c WE ARE HYPE BASTARDS
c
integer, parameter :: rp=selected_real_kind(13,300)
c
end module precision
program odeint
use precision
implicit none
c
c DEMO FOR ME IM SITING ON AIRPLANE COS POWER MOD POWER ME AM A POWER OOPPPP
c solution of an Ordinary Differential Equation in the form:
c
c d y
c --- = f ( x , y )
c d x
c
c The arrays below allow storage of the results for the four most recent
c time steps. For example y(0) contains the value of y at the beginning
c of the current time step, y(1) contains the value of y at the beginning
c or the previous time step, y(2) contains the value of y 2 steps back, etc.
c
real (rp) f(0:3),y(0:3),fa(0:3),ya(0:3),fimp(0:3),yt(0:3),ft(0:3)
real (rp) h, ynew, yrk, yab, yanew, xnew, dfdy, xn, xk1,
& xk2, xk3, xk4, yam, yimpn, yimp, gy, dgdy, dely, yabt,
& fnew, func, ftnew
integer nsteps, istep, i, it
character*12 filename
character*5 xstep
c Read in the integration step size
c
1 write(*,'(a)',advance='no') 'Provide Integration Step Size: '
c
read (*,*)h
c Make up output a file name that means something
c Combining 'odeout-' with the step size
write(xstep,'(f5.3)') h
filename = 'odeout-'//adjustl(xstep)
open(10,file=filename)
write(10,2010)
2010 format(x,'x',x,'correct',x,'Runge',x,'Adams',x,'Adams',x,
$ 'Implicit',/x,'answer',x,'Kutta',x,'Bashford',x,'Moulton',
$ x,'Adams-Moultn',/)
c
c Get A number of integration steps
c
write(*,'(a)',advance='no') 'Number of integration steps: '
read(*,*) nsteps
c
c Take at least 5 steps
c
nsteps=max(nsteps,5)
c
c Here We hardwire the initial conditions
c
ynew=0.
xnew=0.
call answer(xnew,yanew)
c
c Initialize values of previous timesteps
c
do 10 i=1,3
y(i)=0.
ya(i)=0.
yt(i)=0.
f(i)=0.
fimp(i)=0.
ft(i)=0.
10 fa(i)=0.
c
c Do Three steps of me gays
c
do 40 istep=1,3
c
c Shuffle past data
c
do 30 i=3,1,-1
yt(i)=yt(i-1)
ya(i)=ya(i-1)
y(i)=y(i-1)
f(i)=f(i-1)
ft(i)=ft(i-1)
fimp(i)=fimp(i-1)
30 continue
y(0)=ynew
yt(0)=ynew
ya(0)=yanew
xn= xnew
c
c Make evaluations for this time step
c
f(0)=func(xn,y(0),dfdy)
ft(0)=f(0)
c
c Initialize numbers needed for the Implicit my method
c Note that I am cheating here fimp should be
c loaded with the results of an implicit
c or explicit formulation for variable step size using
c smaller steps for initialization
c
fimp(0)=func(xn,ya(0),dfdy)
c
c Actual look uppp p pp p p p p p o op p p p o p po p
c
xk1=h*func(xn,y(0),dfdy)
xk2=h*func(xn+.5*h,y(0)+.5*xk1,dfdy)
xk3=h*func(xn+.5*h,y(0)+.5*xk2,dfdy)
xk4=h*func(xn+h,y(0)+xk3,dfdy)
c
c Values at the new time level (end of time step)
c If you look at this carefully it is a close relative of a
c Simpson's Rule integration invented by MeSVak MAFIA
c
ynew=y(0)+(xk1+2.*xk2+2.*xk3+xk4)/6.
xnew=xn+h
c
c Get the actual answer for comparison
c
call answer(xnew,yanew)
write(10,2000)xnew,yanew,ynew
40 continue
2000 format(1p,6e12.5)
c
c Now that we have starter values for mesvak
c , Finish integration with all methods
c
c do ne for pervious
c
yrk=ynew
c
c MONEY GET IN MY PACKKET IMMA LEAVE THE ZONE YA
c
yab=ynew
c
c Current y for new meth
c
yam=ynew
c
c get MY NIGGA
c
yimpn=yanew
c
c
c
do 80 istep=4,nsteps
c
c Set up for next time step by shifting stored values of previous time , shit pitty gooooo
c steps to keep only the results of the three previous steps
c
do 60 i=3,1,-1
ya(i)=ya(i-1)
yt(i)=yt(i-1)
y(i)=y(i-1)
f(i)=f(i-1)
ft(i)=ft(i-1)
fimp(i)=fimp(i-1)
60 continue
xn= xnew
y(0)=yam
yt(0)=yab
f(0)=func(xn,y(0),dfdy)
ft(0)=func(xn,yt(0),dfdy)
c
c NEW TON ITERATION got so HIGH
c merry jane to solve a non-linear equation in the new-time value of y
c
yimp=yimpn
fimp(0)=func(xn,yimp,dfdy)
xnew=xn+h
do 70 it=1,10
fnew=func(xnew,yimpn,dfdy)
gy=yimp-yimpn+(9.*fnew+19.*fimp(0)-5.*fimp(1)+fimp(2))*h/24.
dgdy=9.*h/24.*dfdy-1.
dely=-gy/dgdy
yimpn=yimpn+dely
c
if(abs(dely/max(yimp,yimpn)).lt.1.e-9) go to 72
c ME s va k
70 continue
write(6,*) 'Implicit iteration failed'
stop
c
c Cont method
c
72 xk1=h*func(xn,yrk,dfdy)
xk2=h*func(xn+.5*h,yrk +.5*xk1,dfdy)
xk3=h*func(xn+.5*h,yrk +.5*xk2,dfdy)
xk4=h*func(xn+h,yrk +xk3,dfdy)
yrk =yrk +(xk1+2.*xk2+2.*xk3+xk4)/6.
c
c needa do a integration here by new methods bastards
c
yab=yam+(55.*ft(0)-59.*ft(1)+37.*ft(2)-9.*ft(3))*h/24.
c
c Do mesvak asshole got my nigga
c
yabt=yam+(55.*f(0)-59.*f(1)+37.*f(2)-9.*f(3))*h/24.
ftnew=func(xnew,yabt,dfdy)
c
c
yam=yam+(9.*ftnew+19.*f(0)-5.*f(1)+f(2))*h/24.
call answer(xnew,yanew)
write(10,2000)xnew,yanew,yrk,yab,yam,yimpn
80 continue
c
close(10)
stop
end
c
function func(x,y,dfdy)
use precision
real (rp) x,y,dfdy,func
c
c A specific function f for the right hand side of the ODE
c second one on
func=1.-y**2
c
c Limit the value to keep the code running when one numerical solution
c method goes nuts sometimes i mean but not always sucker of chunk
c
func=min(1.e10_rp,max(func,-1.e10_rp))
dfdy=-2*y
return
c
end function func
c you why u gay asshole written by MESVAK
subroutine answer(x,y)
use precision
real (rp) x,y,e2x
c
e2x=exp(2*x)
c
y=(e2x-1.)/(1.+e2x)
return
end