Skip to content

Commit

Permalink
Nonlinear index 1 DAE work precision and minor investigations with Ra…
Browse files Browse the repository at this point in the history
…dau methods.
  • Loading branch information
JonasBreuling committed Aug 14, 2024
1 parent 66cd760 commit 2bcb613
Show file tree
Hide file tree
Showing 17 changed files with 310 additions and 26 deletions.
17 changes: 17 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -240,6 +240,23 @@ Starting at $t_0 = \sqrt{1 / 2}$, this problem is solved for $atol = rtol = 10^{

![Weissinger_work_precision](https://raw.githubusercontent.com/JonasBreuling/scipy_dae/main/data/img/Weissinger_work_precision.png)

### Nonlinear index 1 DAE - Kvaernø

In a final example, an nonlinear index 1 DAE is investigated as proposed by [Kvaernø](https://doi.org/10.2307/2008502). The system is given by

$$
\begin{aligned}
(\sin^2(\dot{y}_1) + \sin^2(y_2)) (\dot{y}_2)^2 - (t - 6)^2 (t - 2)^2 y_1 e^{-t} &= 0 \\
(4 - t) (y_2 + y_1)^3 - 64 t^2 e^{-t} y_1 y_2 &= 0 ,
\end{aligned}
$$

It has the analytical solution $y_1(t) = t^4 e^{-t}$, $y_2(t) = t^3 e^{-t} (4 - t)$ and $\dot{y}_1(t) = (4 t^3 - t^4) e^{-t}$, $y_2(t) = (-t^3 + (4 - t) 3 t^2 - (4 - t) t^3) e^{-t}$.

Starting at $t_0 = 0.5$, this problem is solved for $atol = rtol = 10^{-(4 + m / 4)}$, where $m = 0, \dots, 32$. The resulting error at $t_1 = 1$ is compared with the elapsed time of the used solvers in the figure below.

![Kvaerno_work_precision](https://raw.githubusercontent.com/JonasBreuling/scipy_dae/main/data/img/Kvaerno_work_precision.png)

## Install

An editable developer mode can be installed via
Expand Down
Binary file modified data/img/Arevalo_work_precision.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified data/img/Brenan_work_precision.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added data/img/Kvaerno_work_precision.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified data/img/Robertson_work_precision.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified data/img/Weissinger_work_precision.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified data/img/knife_edge_work_precision.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
2 changes: 1 addition & 1 deletion examples/daes/one_element_timoshenko_rod.py
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,7 @@ def f(t, vy):
print(f"F0: {F0}")

# solver options
atol = rtol = 1e-3
atol = rtol = 1e-5

##############
# dae solution
Expand Down
2 changes: 1 addition & 1 deletion scipy_dae/integrate/_dae/benchmarks/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,6 @@ def benchmark(t0, t1, y0, yp0, F, rtols, atols, h0s, name, y_ref=None, y_idx=Non
ax.set_xlabel("||y_ref(t1) - y(t1)||")
ax.set_ylabel("elapsed time [s]")

plt.savefig(f"{name}.png", dpi=300)
plt.savefig(f"data/img/{name}_work_precision.png", dpi=300)

plt.show()
Original file line number Diff line number Diff line change
Expand Up @@ -100,4 +100,4 @@ def sol_true(t):
# reference solution
y_ref = sol_true(t1)[0]

benchmark(t0, t1, vy0, vyp0, F, rtols, atols, h0s, "Knife edge", y_ref)
benchmark(t0, t1, vy0, vyp0, F, rtols, atols, h0s, "knife_edge", y_ref)
53 changes: 53 additions & 0 deletions scipy_dae/integrate/_dae/benchmarks/kvaerno/kvaerno.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
import numpy as np
from scipy_dae.integrate._dae.benchmarks.common import benchmark


"""Nonlinear index 1 DAE, see Kvaerno1990.
References:
-----------
Kvaerno1990: https://doi.org/10.2307/2008502
"""
def F(t, y, yp):
y1, y2 = y
yp1, yp2 = yp
return np.array([
(np.sin(yp1)**2 + np.cos(y2)**2) * yp2**2 - (t - 6)**2 * (t - 2)**2 * y1 * np.exp(-t),
(4 - t) * (y2 + y1)**3 - 64 * t**2 * np.exp(-t) * y1 * y2,
])


def true_sol(t):
return (
np.array([
t**4 * np.exp(-t),
(4 - t) * t**3 * np.exp(-t),
]),
np.array([
(4 * t**3 - t**4) * np.exp(-t),
(-t**3 + (4 - t) * 3 * t**2 - (4 - t) * t**3) * np.exp(-t)
])
)


if __name__ == "__main__":
# exponents
m_max = 32
ms = np.arange(m_max + 1)

# tolerances and initial step size
rtols = 10**(-(4 + ms / 4))
atols = rtols
h0s = 1e-3 * np.ones_like(rtols)

# time span
t0 = 0.5
t1 = 1

# initial conditions
y0, yp0 = true_sol(t0)

# reference solution
y_ref = true_sol(t1)[0]

benchmark(t0, t1, y0, yp0, F, rtols, atols, h0s, "Kvaerno", y_ref)
Loading

0 comments on commit 2bcb613

Please sign in to comment.