381 - A16 - Differential Equations - Runge-Kutta

370 days ago by Professor381

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)