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 = 2
if example == 1:
f = 2*y
t0 = 0
y0 = 2
# for higher order methods....the derivative of f with respect to t treating any y terms with implicit differentiation and noting that y' = f
fp = 2*f # actually fp = 2y' = 2f = 4y
fpp = 4*f # actually fpp = 4y' = 4f = 8y
exact = y0*exp(2*t)
end
if example == 2:
f = t - y
t0 = 0
y0 = 2
fp = 1 - f # actually fp = 1 - y' = 1 - f = 1 - (t-y) = 1 - t + y
fpp = -1 + f # actually fpp = -1 + y'
exact = ((t-1)*exp(t)+3)*exp(-t)
end
# set up discrete time steps for each method below
steps = 8
delt = 1/8
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
|
|
#
var('t')
y = function('y')(t)
exact = desolve(diff(y,t) == t-y, 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 doing 3 explicit steps
t1 = t0 + delt
y1 = y0 + f(t=t0,y=y0)*delt
pretty_print((t1,y1))
t2 = t1 + delt
y2 = y1 + f(t=t1,y=y1)*delt
pretty_print((t2,y2))
t3 = t2 + delt
y3 = y2 + f(t=t2,y=y2)*delt
pretty_print((t3,y3))
Traceback (click to the left of this block for traceback) ... NameError: name 't0' is not defined Traceback (most recent call last): pretty_print((t2,y2))
File "", line 1, in <module>
File "/tmp/tmpyBlgKa/___code___.py", line 3, in <module>
t1 = t0 + delt
NameError: name 't0' is not defined
|
# Euler's method using a loop.
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)
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 Taylors method
ys = [y0]
tay2_pts = [(t0,y0)]
for k in range(steps):
tempy = ys[k] + f(t=ts[k],y=ys[k])*delt + fp(t=ts[k],y=ys[k])*delt^2/2
ys.append(tempy)
tay2_pts.append((ts[k+1],tempy))
end
G = points(tay2_pts,color='magenta')
G += plot(exact,(t,t0,t0+steps*delt),color='green')
S = plot(spline(tay2_pts),(t,t0,t0+steps*delt),color='blue')
show(G+S)
|
# Third order Taylors method
ys = [y0]
tay2_pts = [(t0,y0)]
for k in range(steps):
tempy = ys[k] + f(t=ts[k],y=ys[k])*delt + fp(t=ts[k],y=ys[k])*delt^2/2 + fpp(t=ts[k],y=ys[k])*delt^3/6
ys.append(tempy)
tay2_pts.append((ts[k+1],tempy))
end
G = points(tay2_pts,color='magenta')
G += plot(exact,(t,t0,t0+steps*delt),color='green')
S = plot(spline(tay2_pts),(t,t0,t0+steps*delt),color='blue')
show(G+S)
|
|
|