3

I'm trying to use Python to numerically solve a system of equations described in this paper, Eqs. 30 and 31, with a simplified form looking like:

enter image description here

where G(k) and D(k) are some known functions, independent of Y. Of course, all quantities are functions of t as well. The authors comment that, due to the dependence exhibited by the various functions, a numerical solution is necessary.

I usually implement the solution to this type of coupled equations as indicated here or here, for example, but now the extra k-dependence is confusing me a bit.

Any suggestions? Thanks a lot.

surrutiaquir
  • 93
  • 7
  • 1
    What have you tried, what issue have you encountered? Tell us more about this k coefficient and how it prevent you to solve this problem? – jlandercy Aug 18 '21 at 21:16
  • @jlandercy k is not a parameter but another variable; I have edited the question to show `S` explicitly. Honestly, I don't even know how to start. – surrutiaquir Aug 18 '21 at 21:29
  • Do you need to calculate `f(t)` and `Y(t)` at discrete values of `k` or over a given range? – Patol75 Aug 19 '21 at 01:06
  • @Patol75 ‘k’ is a continuous variable and it will be integrated over a certain range. – surrutiaquir Aug 19 '21 at 02:35
  • I assume we know `D(k)`, right? Does it depend on `Y`? – Patol75 Aug 19 '21 at 02:53
  • @Patol `D(k)` is know and it is Y-independent. – surrutiaquir Aug 19 '21 at 02:59
  • I think what I would do is discretise the range over which `k` is defined, say 100 values. Then, for each value of `k`, apply `solve_ivp` to the first equation. This will give you `f` continuous in `t` but discontinuous in `k`. Then, apply `solve_ivp` to the second equation and use Simpson formula (implemented in SciPy) to obtain `S`. That should work. – Patol75 Aug 19 '21 at 03:00

1 Answers1

1

IDESolver is a general-purpose numerical integro-differential equation solver created by Josh Karpel. Its latest version allows the user to solve multidimensional, coupled IDEs. From the examples provided, an IDE like

enter image description here

with analytical solution (sin x, cos x), can be solved using the following piece of code:

import numpy as np
import matplotlib.pyplot as plt
from idesolver import IDESolver

solver = IDESolver(
            x=np.linspace(0, 7, 100),
            y_0=[0, 1],
            c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]],
            d=lambda x: -0.5,
            f=lambda y: y,
            lower_bound=lambda x: 0,
            upper_bound=lambda x: x,
)

solver.solve()

fig = plt.figure(dpi = 600)
ax = fig.add_subplot(111)

exact = [np.sin(solver.x), np.cos(solver.x)]

ax.plot(solver.x, solver.y[1], label = 'IDESolver Solution', linestyle = '-', linewidth = 3)
ax.plot(solver.x, exact[1], label = 'Analytic Solution', linestyle = ':', linewidth = 3)

ax.legend(loc = 'best')
ax.grid(True)

ax.set_title(f'Solution for Global Error Tolerance = {solver.global_error_tolerance}')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y(x)$')

plt.show()
surrutiaquir
  • 93
  • 7