-
Notifications
You must be signed in to change notification settings - Fork 2
/
jackiewicz.py
64 lines (51 loc) · 1.31 KB
/
jackiewicz.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
55
56
57
58
59
60
61
62
63
64
import time
import numpy as np
import matplotlib.pyplot as plt
from scipy_dae.integrate import solve_dae
"""Jackiewicz's implicit differential equation, see Jackiewicz1981.
References:
-----------
Jackiewicz1981: https://eudml.org/doc/15186
"""
def F(t, y, yp):
return yp - (
(np.sin(t**2 * yp)) / 16
- np.sin(np.exp(y)) / 16
+ 1 / t
)
def true_sol(t):
return np.atleast_1d(np.log(t)), np.atleast_1d(1 / t)
if __name__ == "__main__":
# time span
t0 = 1
t1 = 4
t_span = (t0, t1)
# method = "BDF"
method = "Radau"
# initial conditions
y0, yp0 = true_sol(t0)
# solver options
atol = rtol = 1e-6
# run the solver
start = time.time()
sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method)
end = time.time()
t = sol.t
y = sol.y
success = sol.success
status = sol.status
message = sol.message
print(f"message: {message}")
print(f"elapsed time: {end - start}")
print(f"nfev: {sol.nfev}")
print(f"njev: {sol.njev}")
print(f"nlu: {sol.nlu}")
# visualization
fig, ax = plt.subplots()
ax.set_xlabel("t")
ax.set_ylabel("y")
ax.plot(t, true_sol(t)[0], "-ok", label="y_true")
ax.plot(t, y[0], "--xr", label=f"y {method}")
ax.grid()
ax.legend()
plt.show()