Problem with numerical solving 2nd order nonlinear differential equation

Question

I have this 2nd order nonlinear ODE: $$ \frac{I_0}{C} \left[ e^{\frac{1}{2nU_T} \left(U_{in}(t) - U_{out}(t) - RCU'_{out}(t) - \frac{R}{R_z} U_{out}(t) -LCU''_{out}(t) - \frac{L}{R_z} U'_{out}(t) \right)} - 1 \right] -\frac{U_{out}(t)}{CR_z} = U'_{out}(t). $$

Everything is known except function $U_{out}(t)$. To solve that ODE, I thought about Euler method. So I express $U''_{out}(t)$ as $$ U''_{out}(t) \approx \frac{U'_{out}(t+dt) - U'_{out}(t)}{dt}, $$ then plugin that in equation and express $U'_{out}(t+dt)$ to get $$ U'_{out}(t+dt) = - \frac{2\,n\,dt\,U_T}{LC} \ln \left( \frac{U_{out}(t)}{I_0 R_z} + \frac{C}{I_0}U'_{out}(t) + 1 \right) + (\text{etc, not important}) $$ Then if I express $U_{out}(t+dt)$ as $ U_{out}(t+dt) = U'_{out}(t)dt + U_{out}(t) $, I can with IC $U_{out}(0) = 0, U'_{out}(0)=0$ write a python program to calculate solution. But problem is after some iterations, the argument of $\ln$ becomes negative and program (of course) fails. So, the question is, if there is another method to solve that ODE? thank you

Answer 1

Given the Cauchy problem:

$$ \begin{cases} \frac{I_0}{C}\left[\exp\left(\frac{f(t) - x(t) - RC\dot{x}(t) - \frac{R}{R_z} x(t) -LC\ddot{x}(t) - \frac{L}{R_z} \dot{x}(t)}{2nU_T}\right) - 1\right] -\frac{x(t)}{CR_z} = \dot{x}(t) \\ x(t_0) = x_0 \\ \dot{x}(t_0) = \dot{x}_0 \end{cases} $$

with $t \in [t_0,t_1]$, first of all I would lead to a system of first order ODEs:

$$ \begin{cases} \dot{x}(t) = y(t) \\ \frac{I_0}{C}\left[\exp\left(\frac{f(t) - x(t) - RCy(t) - \frac{R}{R_z} x(t) -LC\dot{y}(t) - \frac{L}{R_z} y(t)}{2nU_T}\right) - 1\right] -\frac{x(t)}{CR_z} = y(t) \\ x(t_0) = x_0 \\ y(t_0) = y_0 \end{cases} $$

so I would write it in normal form:

$$ \begin{cases} \dot{x}(t) = y(t) \\ \dot{y}(t) = \frac{f(t)-2nU_T\ln\left[\frac{C}{I_0}\left(y(t)+\frac{x(t)}{CR_z}\right)+1\right]-x(t)-RCy(t)-\frac{R}{R_z}x(t)-\frac{L}{R_z}y(t)}{LC} \\ x(t_0) = x_0 \\ y(t_0) = y_0 \end{cases}\,. $$

Therefore, having decided the number of intervals $N$, we can apply Euler's method:

$$ \begin{cases} t_k = t_{k-1} + \frac{t_1-t_0}{N} \\ x_k = x_{k-1} + y_{k-1}\frac{t_1-t_0}{N} \\ y_k = y_{k-1} + \frac{f(t_{k-1})-2nU_T\ln\left[\frac{C}{I_0}\left(y_{k-1}+\frac{x_{k-1}}{CR_z}\right)+1\right]-x_{k-1}-RCy_{k-1}-\frac{R}{R_z}x_{k-1}-\frac{L}{R_z}y_{k-1}}{LC}\frac{t_1-t_0}{N} \\ \end{cases} $$

with $k = 1,2,\dots, N$.

All that remains is to draw the graph of $x(t)$ by plotting the points of coordinates $(t_k,x_k)$.

---

Note - The starting ODE:

$$ \frac{I_0}{C}\left[\exp\left(\frac{f(t) - x(t) - RCy(t) - \frac{R}{R_z} x(t) -LC\dot{y}(t) - \frac{L}{R_z} y(t)}{2nU_T}\right) - 1\right] -\frac{x(t)}{CR_z} = y(t) $$

is equivalent to:

$$ \exp\left(\frac{f(t) - x(t) - RCy(t) - \frac{R}{R_z} x(t) -LC\dot{y}(t) - \frac{L}{R_z} y(t)}{2nU_T}\right) = \frac{C}{I_0}\left(y(t) + \frac{x(t)}{CR_z}\right) + 1 $$

which is verified provided that:

$$ \frac{C}{I_0}\left(y(t) + \frac{x(t)}{CR_z}\right) + 1 > 0 $$

and exactly under this condition it's equivalent to writing:

$$ \frac{f(t) - x(t) - RCy(t) - \frac{R}{R_z} x(t) -LC\dot{y}(t) - \frac{L}{R_z} y(t)}{2nU_T} = \ln\left[\frac{C}{I_0}\left(y(t) + \frac{x(t)}{CR_z}\right) + 1\right] $$

i.e.

$$ \dot{y}(t) = \frac{f(t)-2nU_T\ln\left[\frac{C}{I_0}\left(y(t)+\frac{x(t)}{CR_z}\right)+1\right]-x(t)-RCy(t)-\frac{R}{R_z}x(t)-\frac{L}{R_z}y(t)}{LC}. $$

---

Numerical example - Writing in Mathematica:

{c, i0, l, n, r, rz, ut} = {1, 1, 1, 1, 1, 1, 1};
{NN, t0, t1, x0, y0} = {100, 2., 10, 2, -2};
f[t_] := t Cos[t]^2

Δt = (t1 - t0) / NN;
t = x = y = ConstantArray[0, NN + 1];
{t[[1]], x[[1]], y[[1]]} = {t0, x0, y0};

Do[t[[k]] = t[[k - 1]] + Δt;
   x[[k]] = x[[k - 1]] + y[[k - 1]] Δt;
   y[[k]] = y[[k - 1]] + (f[t[[k - 1]]] - 2 n ut Log[c/i0 (y[[k - 1]] + x[[k - 1]]/(c rz)) + 1] 
            - x[[k - 1]] - r c y[[k - 1]] - r/rz x[[k - 1]] - l/rz y[[k - 1]])/(l c) Δt,
  {k, 2, NN + 1}];

ListPlot[Transpose[{t, x}], AxesLabel -> {"t", "x(t)"}, AxesOrigin -> {0, 0}]

we get:

On the other hand, setting y0=-3 doesn't get anywhere and this regardless of the method applied.

---

Addendum - Note the OP's numerical values, I begin the investigation with Mathematica:

{c, i0, l, n, r, rz, ut} = {0.001, 0.0001, 2, 1, 5, 10, 0.0259};
{t0, t1, x0, x0p} = {0., 3, 0, 0};
f[t_] := Sin[2 π t]

xsol = NDSolveValue[{i0/c (Exp[(f[t] - x[t] - r c x'[t] - r/rz x[t] - l c x''[t] - l/rz x'[t])/
                     (2 n ut)] - 1) - x[t]/(c rz) == x'[t], x[t0] == x0, x'[t0] == x0p}, x, {t, t0, t1},
                    Method -> {"EquationSimplification" -> "Residual"}];

Plot[xsol[t], {t, t0, t1}, AxesLabel -> {"t", "x(t)"}, AxesOrigin -> {0, 0}]

So, due to the intrinsic numerical difficulties, I opted for the finite difference method:

$$ \small \begin{cases} \frac{I_0}{C}\left[\exp\left(\frac{f(t) - x(t) - RC\frac{x(t)-x(t-\Delta t)}{\Delta t} - \frac{R}{R_z}x(t) -LC\frac{x(t)-2x(t-\Delta t)+x(t-2\Delta t)}{(\Delta t)^2} - \frac{L}{R_z}\frac{x(t)-x(t-\Delta t)}{\Delta t}}{2nU_T}\right) - 1\right] -\frac{x(t)}{CR_z} = \frac{x(t)-x(t-\Delta t)}{\Delta t} \\ x(t_0) = x_0 \\ \dot{x}(t_0) = \dot{x}_0 \end{cases} $$

that is, by setting:

$$ \begin{aligned} & \alpha \equiv \frac{1}{2nU_T}\left(f(t)+\frac{L\Delta t+CR_z(2L+R\Delta t)}{R_z(\Delta t)^2}x(t-\Delta t)-\frac{CL}{(\Delta t)^2}x(t-2\Delta t)\right); \\ & \beta \equiv \frac{1}{2nU_T}\left(1+\frac{(CR_z+\Delta t)(L+R\Delta t)}{R_z(\Delta t)^2}\right); \\ & \gamma \equiv 1 - \frac{C}{I_0\Delta t}x(t-\Delta t)\,; \\ & \delta \equiv \frac{CR_z+\Delta t}{I_0R_z\Delta t}\,; \\ \end{aligned} $$

it's about solving the equation:

$$ \exp\left[\alpha-\beta x(t)\right] = \gamma+\delta x(t) \quad \Leftrightarrow \quad x(t) = -\frac{\gamma}{\delta}+\frac{1}{\beta}\mathcal{W}\left[\frac{\beta}{\delta}\exp\left(\alpha+\frac{\beta}{\delta}\gamma\right)\right] $$

i.e. it's necessary to iterate:

$$ \begin{cases} t_k = t_{k-1} + \Delta t \\ x_k = -\frac{\gamma}{\delta}+\frac{1}{\beta}\mathcal{W}\left[\frac{\beta}{\delta}\exp\left(\alpha+\frac{\beta}{\delta}\gamma\right)\right] \end{cases} $$

with $\Delta t \equiv \frac{t_1-t_0}{N}$ and $k=2,3,\dots,N$.

For example, in Mathematica we have:

{c, i0, l, n, r, rz, ut} = {0.001, 0.0001, 2, 1, 5, 10, 0.0259};
{NN, t0, t1, x0, x0p} = {1000, 0., 3, 0, 0};
f[t_] := Sin[2 π t]

Δt = (t1 - t0) / NN;
t = x = ConstantArray[0, NN + 1];
{t[[1]], t[[2]]} = {t0, t0 + Δt};
{x[[1]], x[[2]]} = {x0, x0 + x0p Δt};

Do[t[[k]] = t[[k - 1]] + Δt;
   α = 1/(2 n ut) (f[t[[k]]] + (l Δt + c rz (2 l + r Δt))/(rz Δt^2) 
       x[[k - 1]] - c l/Δt^2 x[[k - 2]]);
   β = 1/(2 n ut) (1 + (c rz + Δt) (l + r Δt)/(rz Δt^2));
   γ = 1 - c/(i0 Δt) x[[k - 1]];
   δ = (c rz + Δt)/(i0 rz Δt);
   x[[k]] = -γ/δ + 1/β ProductLog[β/δ Exp[α + β/δ γ]], 
  {k, 3, NN + 1}];

ListPlot[Transpose[{t, x}], AxesLabel -> {"t", "x(t)"}, AxesOrigin -> {0, 0}]

Answer 2

In connection with the response from @TeM, here are numerical values for parameters:

\begin{align} U_{in}(t) &= U_{in_A} \sin 2 \pi f t \\ U_{in_A} &= 1 \\ f &= 1 \\ C &= 0.001 \\ R &= 5 \\ R_{Z} &= 10 \\ L &= 2 \\ I_{0} &= 0.0001 \\ n &= 1 \\ U_T &= 0.0259 \\ \end{align} With these values the program fails (With values from @TeM program works). I can change the parameters $U_{in_A}, f, C, R, R_Z, L$. Values of $I_0, U_T$ must be of the same order (can't be much bigger or smaller than now)

Answer 3

Too messy for a comment

With the following MATHEMATICA script, for trivial initial conditions $U_{out}(0) = U'_{out}(0)=0$ we obtain some results.

Ui[t_] := Ui0 Sin[2 Pi f t]
tmax = 3;
parms = {Ui0 -> 1, f -> 1, C0 -> 0.001, R -> 5, Rz -> 10, L -> 2, I0 -> 0.0001, n -> 1, UT -> 0.0259};
ode = I0/C0 (Exp[1/(2 n UT) (Ui[t] - Uo[t] - R C0 Uo'[t] - R/Rz Uo[t] - 
      L C0 Uo''[t] - L/Rz Uo'[t])] - 1) - Uo[t]/C0/Rz - Uo'[t] /. parms
solUo = NDSolve[{ode == 0, Uo[0] == Uo'[0] == 0}, Uo, {t, 0, tmax}][[1]]

Uot = Uo[t] /. solUo;
Plot[Uot, {t, 0, 3}, PlotStyle -> {Thick, Blue}]

So perhaps using a higher order differential operator, the python procedure could work.

NOTE

MATHEMATICA handles this ODE as a differential-algebraic system

$$ \cases{K_0\left(e^{y_1(t)}-1\right)= y_2(t)\\ \mathcal{D}_1U_{out}(t) = y_1(t)\\ \mathcal{D}_2U_{out}(t) = y_2(t) } $$

where $\mathcal{D}_1[],\mathcal{D}_2[]$ are linear differential operators.