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