Heun's Method p. 110 #4

818 days ago by dwatson

Draw a directional field to get an idea of what the solution is:

t,y=var('t,y') df = 3*cos(t)-2*y v=plot_slope_field(df,(t,0,3),(y,-3,3)) plot(v) 
       

#4a. Find approximate values of the solution of the given initial value problem at t = 0.1, 0.2, 0.3, and 0.4 using the Euler method with h = 0.1

# differential equation f(t,y) = 3*cos(t)-2*y g(t,y) = diff(f(t,y),t) # initial value t = 0 y = 0 # step size h = 0.1 # Heun's method y_n+1 = y_n + 1/2*(f_n+f_n+1)*h for j in range(4) : ybar = y+f(t,y)*h #approximation of y_n+1 y = y + 1/2*(f(t,y)+f(t+h,ybar))*h t = t + h print t.n(digits=2), y.n(digits=4) 
       
0.10 0.2692
0.20 0.4872
0.30 0.6604
0.40 0.7943
0.10 0.2692
0.20 0.4872
0.30 0.6604
0.40 0.7943

#4b. Repeat part (a) with h = 0.05.  Compare the results with those found in (a).

f(t,y) = 3*cos(t)-2*y g(t,y) = diff(f(t,y),t) t = 0 y = 0 h = 0.05 for i in range(4) : for j in range(2) : ybar = y+f(t,y)*h y = y + 1/2*(f(t,y)+f(t+h,ybar))*h t = t + h print t.n(digits=2), y.n(digits=5) 
       
0.10 0.27092
0.20 0.49003
0.30 0.66403
0.40 0.79847
0.10 0.27092
0.20 0.49003
0.30 0.66403
0.40 0.79847

#4c. Repeat part (a) with h = 0.025.  Compare the results with those found in (a) and (b).

f(t,y) = 3*cos(t)-2*y t = 0 y = 0 h = 0.025 for i in range(4) : for j in range(4) : ybar = y+f(t,y)*h y = y + 1/2*(f(t,y)+f(t+h,ybar))*h t = t + h print t.n(digits=2), y.n(digits=5) 
       
0.10 0.27131
0.20 0.49069
0.30 0.66487
0.40 0.79942
0.10 0.27131
0.20 0.49069
0.30 0.66487
0.40 0.79942

#4d. Find the solution of the given problem and evaluate at t = 0.1, 0.2, 0.3, and 0.4.  Compare these values with the results of (a), (b), and (c).

x = var('t') y = function('y', t) a = desolve(diff(y,t)-3*cos(t)+2*y, y, ics=[0,0]) a.expand() 
       
6/5*cos(t) - 6/5*e^(-2*t) + 3/5*sin(t)
6/5*cos(t) - 6/5*e^(-2*t) + 3/5*sin(t)
t = 0.1 h = 0.1 for j in range(4) : a(t) t = t + h 
       
0.271428144628150
0.490897436643760
0.665141947634699
0.799729441247987
0.271428144628150
0.490897436643760
0.665141947634699
0.799729441247987
d = contour_plot(C, (x,-2, 2), (y,-2,2),fill = False, linewidths = 5, contours = 30, cmap = 'hsv')