352 - Topic 04 - Basic Numerical Solutions

# Worksheet still under development. Use at your own risk! 8-o 

reset() 

t,y=var('t y') f(t,y) = (y^2-2) - t # ff(t,z) = z^2-2 # f(t,y) = t^3*y+y^3*t-2*y^2 # ff(t,z) = t^3*z+z^3*t-2*z^2 # f(t,y) = t-y # ff(t,z) = t-z # f(t,y) = y^2-t^2 # ff(t,z) = z^2-t^2 # z = function('z') # soln = function('soln') # soln = desolve(diff(z,t) == ff(t,z), z) # print "The General Solution is" # show(soln) # show(solve(soln,z)) 

print f 

(t, y) |--> y^2 - t - 2

(t, y) |--> y^2 - t - 2

# Euler's Method show("y' = ",f) @interact def _(t0=input_box(0,label='$t_0$',width=10), tend = input_box(1,label='$t_{end}$',width=10), y0=input_box(1,label='$y_0$',width=10), delt=[1/2,1/5,1/10,1/20,1/50], show_soln = checkbox(default=false)): steps = floor((tend-t0)/delt) w = y0 t1 = t0 G = Graphics() G += plot_slope_field(f(t,y),(t,t0,tend),(y,-1,2),color='gray') for k in range(steps): w1 = w + delt*f(t1,w) G += line(((t1,w),(t1+delt,w1)),color='red',thickness=3,title = "Euler's Method") G += point((t1+delt,w1),color='green',size=50,zorder=3) w = w1 t1 = t1+delt # G += desolve_rk4(f(t,y),y,ics=[t0,y0],ivar=t,output='slope_field',end_points=[t0,t0+steps*delt],thickness=3) # Trying to play with the exact symbolic solution # soln = desolve(diff(z,t) - ff(z,t),z,ics=(t0,y0)) # print ff,soln # G += plot(soln,(t,t0,tend)) show(G) 

# Modified Euler's (Heun's) @interact def _(t0=input_box(0,label='$t_0$',width=10), tend = input_box(1,label='$t_{end}$',width=10), y0=input_box(1,label='$y_0$',width=10), delt=[1/2,1/5,1/10,1/20,1/50], show_soln = checkbox(default=false)): delt2 = delt/2 steps = floor((tend-t0)/delt) w = y0 t1 = t0 G = Graphics() for k in range(steps): t2 = t1+delt w1a = w + delt*f(t1,w) w1 = w + delt2*(f(t1,w)+f(t2,w1a)) G += line(((t1,w),(t2,w1)),color='red',thickness=3,title = "Modified Euler's Method") G += point((t1+delt,w1),color='green',size=50,zorder=3) w = w1 t1 = t1+delt if show_soln: G += desolve_rk4(f(t,y),y,ics=[t0,y0],ivar=t,output='slope_field',end_points=[t0,t0+steps*delt],thickness=3) show(G) 

# Runge-Kutta Method @interact def _(t0=input_box(0,label='$t_0$',width=10), tend = input_box(1,label='$t_{end}$',width=10), y0=input_box(1,label='$y_0$',width=10), delt=[1/2,1/5,1/10,1/20,1/50], show_soln = checkbox(default=false)): delt2 = delt/2 steps = floor((tend-t0)/delt) w = y0 t1 = t0 G = Graphics() for k in range(steps): t2 = t1+delt k1 = f(t1,w) k2 = f(t1+delt2,w+delt*k1/2) k3 = f(t1+delt2,w+delt*k2/2) k4 = f(t1+delt,w+delt*k3) w1 = w + (delt/6)*(k1+2*k2+2*k3+k4) G += line(((t1,w),(t2,w1)),color='red',thickness=3,title = "Runge Kutta Method") G += point((t1+delt,w1),color='green',size=50,zorder=3) w = w1 t1 = t1+delt if show_soln: G += desolve_rk4(f(t,y),y,ics=[t0,y0],ivar=t,output='slope_field',end_points=[t0,t0+steps*delt],thickness=3) # Notice that Sage "solves" this DE numerically using this exact 4th order Runge-Kutta method anyway. show(G) 

# still working on this one below too. 

# add a visual that shows several numerical solution curves using progressively # smaller step sizes. @interact def _(t0=input_box(0,label='$t_0$',width=10), y0=input_box(1,label='$y_0$',width=10)): w = y0 t1 = t0 G = Graphics() for delt in (1,1/2,1/4,1/10,1/100): steps = 2/delt for k in range(steps): w1 = w + delt*f(t1,w) G += point((t1+delt,w1),color='green',size=50,zorder=3) G += line(((t1,w),(t1+delt,w1)),color='red',thickness=3) w = w1 t1 = t1+delt # G += desolve_rk4(f(t,y),y,ics=[t0,y0],ivar=t,end_points[t0,t0+2],thickness=3) show(G)