4
$\begingroup$

Hi I would like to solve the following equations:

$\frac{dx(t)}{dt} = \gamma v- k x(t) \left(1-b e^{f(t)}\left(1+b e^{f(t)}\right)\right)$

where

$f(t) = k x(t) \left(1+b e^{f(t)}\right)$

As u see this the solution is self consistent, meaning that for each value of $x$ I need to determine first $f(t)$, plug it inside the equation and then solve.

Any suggestions or reference about how to use Findroot or any other method in NDSolve?

$\endgroup$
4
$\begingroup$

NDSolve can solve differential algebraic equations (DAEs) directly. Following @AccidentalFourierTransform's simplifying assumptions,

γ = 0;
b = k = 1;

sol = NDSolve[{
  x'[t] == γ v - k x[t] (1 - b Exp[f[t]] (1 + b Exp[f[t]])),
  f[t] == k x[t] (1 + b Exp[f[t]]),
  x[0] == -0.5}, {x, f}, {t, 0, 1}][[1]];

Plot[x[t] /. sol, {t, 0, 1}]

Mathematica graphics

In your particular example, you can explicitly solve for f[t] using Solve and then plug that into NDSolve.

Clear[k, b]
fsol = Solve[f[t] == k x[t] (1 + b Exp[f[t]]), f[t]][[1]]
(* {f[t] -> -ProductLog[-b E^(k x[t]) k x[t]] + k x[t]} *)

γ = 0;
b = k = 1;

sol = NDSolve[{
  x'[t] == γ v - k x[t] (1 - b Exp[f[t]] (1 + b Exp[f[t]])),
  x[0] == -0.5} /. fsol, {x}, {t, 0, 1}][[1]];

gives an identical result.

$\endgroup$
4
$\begingroup$

Take $\gamma=0$ and $k=b=1$ to simplify the problem. The general case is then straightforward to implement.

Define

f[x_?NumericQ] := FindRoot[z == x (1 + Exp[z]), {z, 1}][[1, 2]]

With this, and using $x(0)=-0.5$ as the initial condition, we get

NDSolve[{x'[t] == -x[t] (1 - Exp[f[x[t]]] (1 + Exp[f[x[t]]])), x[0] == -.5}, x, {t, 0, 1}]
Plot[%[[1, 1, 2]][t], {t, 0, 1}]

enter image description here

Ta-da!

$\endgroup$
  • $\begingroup$ I think you need minus in your x'[t] equation to match OP's system. $\endgroup$ – Chris K Jan 14 '18 at 17:48
  • $\begingroup$ @ChrisK yup, thanks! $\endgroup$ – AccidentalFourierTransform Jan 14 '18 at 17:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.