1

I'm new to the SciPy.org maths libraries, so this may be a fairly basic question for those familiar with them.

For this ODE:

y'(t) - 0.05y(t) = d, y(0) = 10

how do I calculate the value of 'd' if y(10) = 100?

I can solve for y(t) this way:

import sympy as sym
y = sym.Function('y')
t, d = sym.symbols('t d')
y1 = sym.Derivative(y(t), t)
eqdiff = y1 - 0.05*y(t) - d
sol = sym.dsolve(eqdiff, y(t), ics={y(0): '10'})
sol

y(t)= −20.0d + (20.0d + 10.0)e^(0.05t)

Whether "sol" is usable to solve for d when y(10) = 100 is unknown to me (SymPy may not be the library of choice for this).

I've looked at numerous web pages such as these for ideas but haven't found a way forward:

https://docs.sympy.org/latest/modules/solvers/ode.html

Converting sympy expression to numpy expression before solving with fsolve( )

https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

I'm aware there are graphical ways to address the problem, but I want a numeric result.

Thanks in advance for helpful advice.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
hilbert
  • 13
  • 3
  • This can be solved by hand rather easily... – Julien Oct 03 '21 at 07:27
  • @Julien Certainly true, but the issue isn't to learn to solve this particular equation. It's to learn the general method for problems that aren't so easy. ; ) – hilbert Oct 03 '21 at 15:28

2 Answers2

4

You can substitute the values and use solve:

In [5]: sol.subs(t, 10)
Out[5]: y(10) = 12.9744254140026⋅d + 16.4872127070013

In [6]: sol.subs(t, 10).subs(y(10), 100)
Out[6]: 100 = 12.9744254140026⋅d + 16.4872127070013

In [7]: solve(sol.subs(t, 10).subs(y(10), 100), d)
Out[7]: [6.43672337141557]

https://docs.sympy.org/latest/modules/solvers/solvers.html#sympy.solvers.solvers.solve

Oscar Benjamin
  • 12,649
  • 1
  • 12
  • 14
2

You can also solve it with scipy. The whole task is a boundary value problem with a free parameter, one state dimension plus one free parameter equals two boundary slots. So use solve_bvp (even if it is a scalar problem, the solver treats every state space as vector space)

def eqn(t,y,d): return d+0.05*y

def bc(y0,y10,d): return [ y0[0]-10, y10[0]-100 ]

x_init = [0,10]
y_init = [[10, 100]]
d_init = 1

res = solve_bvp(eqn, bc, x_init, y_init, p=[d_init], tol=1e-5)
print(res.message, f"d={res.p[0]:.10f}")

which gives

The algorithm converged to the desired accuracy. d=6.4367242595
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51