-
Notifications
You must be signed in to change notification settings - Fork 190
/
ForwardEuler_func.py
54 lines (44 loc) · 1.6 KB
/
ForwardEuler_func.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
"""Function implementing the Forward Euler method for scalar ODEs."""
import numpy as np
def ForwardEuler(f, U0, T, n):
"""Solve u'=f(u,t), u(0)=U0, with n steps until t=T."""
t = np.zeros(n+1)
u = np.zeros(n+1) # u[k] is the solution at time t[k]
u[0] = U0
t[0] = 0
dt = T/float(n)
for k in range(n):
t[k+1] = t[k] + dt
u[k+1] = u[k] + dt*f(u[k], t[k])
return u, t
# Problem: u'=u
def f(u, t):
return u
u, t = ForwardEuler(f, U0=1, T=4, n=20)
# Compare numerical solution and exact solution in a plot
from scitools.std import plot, exp
u_exact = exp(t)
plot(t, u, 'r-', t, u_exact, 'b-',
xlabel='t', ylabel='u', legend=('numerical', 'exact'),
title="Solution of the ODE u'=u, u(0)=1", savefig='tmp.pdf')
# More exact verification
def test_ForwardEuler_against_hand_calculations():
"""Verify ForwardEuler against hand calc. for 3 time steps."""
u, t = ForwardEuler(f, U0=1, T=0.2, n=2)
exact = np.array([1, 1.1, 1.21]) # hand calculations
error = np.abs(exact - u).max()
success = error < 1E-14
assert success, '|exact - u| = %g != 0' % error
def test_ForwardEuler_against_linear_solution():
"""Use knowledge of an exact numerical solution for testing."""
def f(u, t):
return 0.2 + (u - u_exact(t))**4
def u_exact(t):
return 0.2*t + 3
u, t = ForwardEuler(f, U0=u_exact(0), T=3, n=5)
u_e = u_exact(t)
error = np.abs(u_e - u).max()
success = error < 1E-14
assert success, '|exact - u| = %g != 0' % error
test_ForwardEuler_against_hand_calculations()
test_ForwardEuler_against_linear_solution()