481 - Euler's Method

24 days ago by Professor481

John Travis

Mississippi College

# number of time steps to display n = 20 # differential equation t,y=var('t,y') f(t,y) = y*(3-t*y) # f(t,y) = y/2 # initial value t0 = 0 y0 = 1 
       

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

t,y=var('t,y') v=plot_slope_field(f,(t,0,8),(y,-3,5)) 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

eul1 = [(t0,y0)] # step size h = 0.1 tj = t0 w = y0 # Euler's method y_n+1 = y_n + f_n*h for j in range(n) : w = w + f(tj,w)*h tj = tj + h print tj.n(digits=2), w.n(digits=4) eul1.append((tj,w)) 
       
0.10 1.300
0.20 1.673
0.30 2.119
0.40 2.620
0.50 3.131
0.60 3.581
0.70 3.886
0.80 3.994
0.90 3.916
1.0 3.711
1.1 3.447
1.2 3.174
1.3 2.917
1.4 2.686
1.5 2.482
1.6 2.302
1.7 2.145
1.8 2.006
1.9 1.884
2.0 1.775
0.10 1.300
0.20 1.673
0.30 2.119
0.40 2.620
0.50 3.131
0.60 3.581
0.70 3.886
0.80 3.994
0.90 3.916
1.0 3.711
1.1 3.447
1.2 3.174
1.3 2.917
1.4 2.686
1.5 2.482
1.6 2.302
1.7 2.145
1.8 2.006
1.9 1.884
2.0 1.775

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

tj = t0 w = y0 eul2 = [(t0,y0)] h = 0.05 for i in range(n) : for j in range(2) : w = w + f(tj,w)*h tj = tj + h print tj.n(digits=2), w.n(digits=5) eul2.append((tj,w)) 
       
0.10 1.3192
0.20 1.7176
0.30 2.1902
0.40 2.7091
0.50 3.2156
0.60 3.6301
0.70 3.8820
0.80 3.9461
0.90 3.8510
1.0 3.6535
1.1 3.4085
1.2 3.1540
1.3 2.9109
1.4 2.6883
1.5 2.4888
1.6 2.3117
1.7 2.1550
1.8 2.0162
1.9 1.8929
2.0 1.7831
0.10 1.3192
0.20 1.7176
0.30 2.1902
0.40 2.7091
0.50 3.2156
0.60 3.6301
0.70 3.8820
0.80 3.9461
0.90 3.8510
1.0 3.6535
1.1 3.4085
1.2 3.1540
1.3 2.9109
1.4 2.6883
1.5 2.4888
1.6 2.3117
1.7 2.1550
1.8 2.0162
1.9 1.8929
2.0 1.7831

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

tj = t0 w = y0 eul3 = [(t0,y0)] h = 0.025 for i in range(n) : for j in range(4) : w = w + f(tj,w)*h tj = tj + h print tj.n(digits=2), w.n(digits=5) eul3.append((tj,w)) 
       
0.10 1.3300
0.20 1.7423
0.30 2.2288
0.40 2.7551
0.50 3.2561
0.60 3.6508
0.70 3.8775
0.80 3.9231
0.90 3.8217
1.0 3.6276
1.1 3.3901
1.2 3.1433
1.3 2.9062
1.4 2.6878
1.5 2.4909
1.6 2.3153
1.7 2.1592
1.8 2.0206
1.9 1.8973
2.0 1.7872
0.10 1.3300
0.20 1.7423
0.30 2.2288
0.40 2.7551
0.50 3.2561
0.60 3.6508
0.70 3.8775
0.80 3.9231
0.90 3.8217
1.0 3.6276
1.1 3.3901
1.2 3.1433
1.3 2.9062
1.4 2.6878
1.5 2.4909
1.6 2.3153
1.7 2.1592
1.8 2.0206
1.9 1.8973
2.0 1.7872

#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).

t = var('t') y = function('y')(t) a = desolve(diff(y,t)-y*(3-y*t), y, ics=[0,1]) # a = desolve(diff(y,t)-y/2, y, ics=[0,1]) a.expand().show() 
       

                                
                            

                                
tstep = 0 h = 0.1 for i in range(n): a(t=tstep).n() tstep = tstep + h 
       
1.00000000000000
1.34164513165038
1.76882747200711
2.26946257207044
2.80204302139172
3.29513225310198
3.66899339220405
3.87169331802918
3.90086678200998
3.79417826173136
3.60306882930609
3.37226213538943
3.13235535282019
2.90096201317796
2.68660263679334
2.49232228669963
2.31821494061443
2.16295851248213
2.02467138141736
1.90135770805648
1.00000000000000
1.34164513165038
1.76882747200711
2.26946257207044
2.80204302139172
3.29513225310198
3.66899339220405
3.87169331802918
3.90086678200998
3.79417826173136
3.60306882930609
3.37226213538943
3.13235535282019
2.90096201317796
2.68660263679334
2.49232228669963
2.31821494061443
2.16295851248213
2.02467138141736
1.90135770805648
d = plot(a,(x,0, 2.5))+points(eul1,color='red',size=20)+points(eul2,color='blue',size=20)+points(eul3,color='green',size=20) t,y=var('t,y') v=plot_slope_field(f,(t,0,2.5),(y,0,4),plot_points=40) plot(v+d)