Solve numerically the initial value problem:
$ y' = f(t,y) $
$ y(a) = y_0 $
So, below we will need be given $f(t,y)$ and the starting point $(a,y_0)$.
var('t,y')
example = 4
if example == 1:
f = 2*y
t0 = 0
y0 = 2
fp = 2*f # actually fp = 2y' = 2f = 4y
exact = y0*exp(2*t)
end
if example == 2:
f = t - y
t0 = 0
y0 = 2
exact = ((t-1)*exp(t)+3)*exp(-t)
end
if example == 3:
f = -y*cos(t)
t0 = 2
y0 = 1/2
exact = 1/2*exp(sin(2))*exp(-sin(t))
end
if example == 4:
f = 1+y/t
t0 = 1
y0 = 1
exact = t*(log(t)+y0)
end
# set up discrete time steps for each method below
steps = 11
delt = 1/4
ts = [t0]
for k in range(steps+1): # create uniform steps to make things easy. Not using any adjustable step size
ts.append(ts[k] + delt)
end
pretty_print(ts)
|
|
#
y = function('y')(t)
# exact = desolve(diff(y,t) == t-y, y, ics=[t0,y0])
exact = desolve(diff(y,t) == -y*cos(t), y, ics=[t0,y0])
show(exact)
# Need to make this automatic using f and not having to explicitly enter f's formula
# exact = desolve(diff(z,x) == f(t=t,y=z), z, ics=[t0,y0])
# show(exact)
|
|
# Euler's method is the first order RK method.
ys = [y0]
euler_pts = [(t0,y0)]
for k in range(steps):
tempy = ys[k] + f(t=ts[k],y=ys[k])*delt
ys.append(tempy.n()) # all we want is the numerical value, not symbolic
euler_pts.append((ts[k+1],tempy))
end
G = points(euler_pts,color='magenta')
G += plot(exact,(t,t0,t0+steps*delt),color='green')
S = plot(spline(euler_pts),(t,t0,t0+steps*delt),color='blue')
show(G+S)
|
# Second order RK method
# explicit midpoint method
ys = [y0]
RK_pts = [(t0,y0)]
for k in range(steps):
midy = ys[k] + f(t=ts[k],y=ys[k])*delt/2
tempy = ys[k] + f(t=ts[k]+delt/2,y=midy)*delt
ys.append(tempy.n())
RK_pts.append((ts[k+1],tempy))
end
G = points(RK_pts,color='magenta')
G += plot(exact,(t,t0,t0+steps*delt),color='green')
S = plot(spline(RK_pts),(t,t0,t0+steps*delt),color='blue')
show(G+S)
|
# Fourth order RK method
# known as RK4 and is the standard that people use for high accuracy
ys = [y0]
RK_pts = [(t0,y0)]
for k in range(steps):
K1 = f(t=ts[k],y=ys[k])
K2 = f(t=ts[k]+delt/2,y=ys[k]+K1*delt/2)
K3 = f(t=ts[k]+delt/2,y=ys[k]+K2*delt/2)
K4 = f(t=ts[k]+delt,y=ys[k]+K3*delt)
tempy = ys[k] + (K1+2*K2+2*K3+K4)*delt/6
ys.append(tempy.n())
RK_pts.append((ts[k+1],tempy))
end
G = points(RK_pts,color='magenta')
G += plot(exact,(t,t0,t0+steps*delt),color='green')
S = plot(spline(RK_pts),(t,t0,t0+steps*delt),color='blue')
show(G+S)
|
|
|