|
=={{header|Python}}==
#######################################################
# RHS of the ODE to be solved
def yp(t,y):
return t*sqrt(y)
▲<lang Python> def RK4(f, dt):
#######################################################
return lambda t, y: (
# The RK4 routine
lambda dy1: (
#
lambda dy2: (
# INPUT ARGUMENTS
lambda dy3: (
# yp function object providing the RHS of the 1st order ODE
lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
# tmin, tmax the domain of the solution
)( dt * f( t + dt , y + dy3 ) )
# y0, t0 initial values
dy3= )( dt * yp f( tnt + dt/2 .0, yny + dy2/2 .0 ) ) ▼
# dt time step over which the solution is sought
dy2= )( dt * yp f( tnt + dt/2 .0, yny + dy1/2 .0 ) ) ▼
#
)( dt * f( t , y ) )
# OUTPUT
# ti the list of points where the equation has been solved
# yi the values of the solution found
########################################################
def RK4(yp,tmin,tmax,y0,t0,dt):
yi=[y0]
ti=[t0]
nmax=int( (tmax-tmin)/dt +1)
for n in range(1,nmax,1):
tn=ti[n-1]
yn=yi[n-1]
dy1=dt*yp( tn, yn)
▲ dy2=dt*yp( tn+dt/2.0, yn+dy1/2.0)
▲ dy3=dt*yp( tn+dt/2.0, yn+dy2/2.0)
dy4=dt*yp( tn+dt, yn+dy3)
yi.append( yn+(1.0/6.0)*(dy1+2.0*dy2+2.0*dy3+dy4) )
ti.append(tn+dt)
return [ti,yi]
def theory(t): return (t**2 + 4)**2 /16
#######################################################
dt = .1
#
[tmin, tmax, dt] = [0.0, 10.0, 0.1]
dy = RK4(lambda t, y: t*sqrt(y), dt)
[t,y]=RK4(yp,tmin,tmax,y0,t0,dt)
for i in range(0,len(t),10):
▲[y0t, t0]y = [10. 0, 01. 0]
print ("y(%2.1f)\t= %4.6f \t error:%4.6g")%(t[i], y[i], abs( y[i]- ((t[i]**2 + 4.0)**2)/16.0 ) ) ▼while t <= 10:
if abs(round(t) - t) < 1e-5:
▲ print ("y(%2.1f)\t= %4.6f \t error: %4.6g" ) % ( t [i], y [i], abs( y [i] - (theory(t [i]**2 + 4.0) **2) /16.0 ) )
t, y = t + dt, y + dy( t, y )
</lang>
{{Out}}
'''Output'''
<pre>y(0.0) = 1.000000 error: 0
<pre>
y(01.0) = 1.000000562500 error: 01.45722e-07
y(12.0) = 13.562500999999 error:1 9.45722e19479e-07
y(23.0) = 310.999999562497 error:9 2.19479e90956e-0706
y(34.0) = 1024.562497999994 error:2 6.90956e23491e-06
y(45.0) = 2452.999994562489 error:6 1.23491e08197e-0605
y(56.0) = 5299.562489999983 error: 1.08197e65946e-05
y(67.0) = 99175.999983562476 error:1 2.65946e35177e-05
y(78.0) = 175288.562476999968 error:2 3.35177e15652e-05
y(89.0) = 288451.999968562459 error:3 4.15652e07232e-05
y(910.0) = 451675.562459999949 error:4 5.07232e09833e-05</pre>
y(10.0) = 675.999949 error:5.09833e-05
</pre>
=={{header|Racket}}==
|