# Solving coupled first-order autonomous differential equations

I am trying to run a simulation of a circuit to observe chaotic behavior, and have obtained a set of differential equations from an analysis of the circuit.

I thought perhaps I could produce a plot of the current j[t] or generate a parametric plot of j[t] and the voltage V[t]. My constants and used functions I used were

A=1;
ϕ=0.0257275275;
ω=151.6*10^3;
Cj=17*10^(-9);
Cs=35*10^(-9);
Rt=223.9;
I0=5*10^-6;
L=15*10^-3;
V0[t_]:=A*Cos[ω*t];
Id[x_]:=I0*(Exp[x/ϕ]-1);
Cd[x_] := Piecewise[
{{Cj/Sqrt[1 - x/ϕ], x < 0},
{Cs*Exp[x/ϕ], x >= 0}}
];


and the differential equations and input I used were,

eqn = {V'[t] == (j[t] - Id[V[t]])/Cd[V[t]], j'[t] == (V0[t] - Rt*j[t] - V[t])/L};
bc = {V[0] == A, j[0] == 0, θ[0] == 0};
pl = NDSolve[{eqn, bc}, j, {t, 0, 1}]
ParametricPlot[Evaluate[j[t] /. pl], {t, 0, 1}, PlotRange -> All]


It would appear however that my method of evaluation is inappropriate as I'm getting an error from NDSolve: The method currently implemented for delay differential equations does not support delays that depend directly on the time variable or dependent variables. (Resolved)

Is there anyway I can resolve this or perhaps another way to solve the set of differential equations?

Edit: Removed unnecessary third diff. eq. θ'[t]=ω. Unfortunately, it seems to be unable to solve an ordinary differential equation. A research paper on a similar circuit said that the equations could be solved numerically by an explicit fourth-order Runge-Kutta algorithm, so perhaps there's a way I can instruct Mathematica to solve it using this?

Edit 2: RemovedCj/Sqrt[1 - x/ϕ] from x≥0 for Cd

• There's an easy closed form for θ[t]; why not help NDSolve[] out and simplify your DEs? Commented Jun 7, 2015 at 10:55
• @J.M. That was an easy fix, thanks. Now it seems to be having an issue solving the autonomous equations. I've added an approach used in the past from a paper I found in Edit. Commented Jun 7, 2015 at 17:31
• Copying and paste 2D input doesn't always work well (e.g. the piecewise function). It's better to copy-paste InputForm and to check by copying the SE code back into Mathematica and testing it. I edited your Q to fix it, but I had to assume what the piecewise was meant to be. Please check it. Commented Jun 7, 2015 at 18:15
• The difficulty is probably the term V[ω*t], which will require the value of V to be known far in advance of the current time t, in order to compute the value for j'[t]; further the difference is not fixed, but varies, which is what the error is complaining about. Are you certain about the V[ω*t] term? Commented Jun 7, 2015 at 18:17
• I've reported the "DiscontinuityProcessing" issue to Wolfram Support. In a couple of days, we'll see if they think it's a bug. -- The techie at WRI called it an "issue", but I'll take bug as read. Commented Jun 10, 2015 at 20:53

After removing the superfluous θ[0] == 0 from bc, the key to the solution is a subtle change in the use of Piecewise in the definition of Cd:

Cd[x_] := Piecewise[{{Cj/Sqrt[1 - x/ϕ], x < 0}}, Cs*Exp[x/ϕ]]

eqn = {V'[t] == (j[t] - Id[V[t]])/Cd[V[t]], j'[t] == (V0[t] - Rt*j[t] - V[t])/L};
bc = {V[0] == A, j[0] == 0 (*, \[Theta][0] == 0*)};
pl = NDSolve[{eqn, bc}, j, {t, 0, 1}]


With the OP's definition of Cd, we get

Power::infy: Infinite expression 1/0 encountered. >>

NDSolve::nlnum: The function value {31.85,ComplexInfinity} is not a list of numbers with dimensions {2} at {t,j[t],V[t],NDSolves$50801[t],NDSolves$50800[t]} = {0.00612563,-0.000795134,0.,1,1}. >>

Here's how I thought about debugging this one: One should look carefully at the data, which luckily has not been Short-ened:

{t, j[t], V[t], NDSolves$50801[t], NDSolves$50800[t]} =
{0.0061256268269912675, -0.0007951340768516366, 0., 1, 1}


It seems interesting that a time t = 0.0061256268269912675, the variable V[t] is exactly 0.. Two more messages report the same thing at succeeding times, t = 0.00621476 and t = 0.00624896. But why should V[t] = 0., be of importance? The original equations do not divide by V[t].

eqn


Ah, we see two things. First, there is a discontinuity at V[t] == 0, so NDSolve will solve for the event where V[t] == 0. Second, we do divide by the Piecewise function, which has a 0 in it. And the error was 1/0, not 1/0.. That must be it. And it was. (But I'm afraid I don't why the 0 leaked out. Evaluating eqn by hand does not yield an error:

eqn /. Thread[{t, j[t], V[t], NDSolves$50801[t], NDSolves$50800[t]} ->
{0.0061256268269912675, -0.0007951340768516366, 0., 1, 1}]
(*
{Derivative[1][V][0.00612563] == -22718.1,
Derivative[1][j][0.00612563] == 31.85}
*)


Maybe it's a bug somewhere in the event handling of the discontinuity processing of NDSolve.

If we solve either version, the OP's or mine,

plOP = NDSolve[{eqn, bc}, {j, V}, {t, 0, 0.01}]
pl = NDSolve[{eqn, bc}, {j, V}, {t, 0, 0.01}]
plOP === pl
(*  True  *)


(sorry but be sure to change the Cd and eqn), we get identical solutions even though the OP's generate 1/0 errors, which are normally catastrophic. We can see that V[t] hit 0. around t == 0.006 and then starts to "chatter." Whether or not it makes physical sense, probably the integration should stop there. (If it does not make physical sense, then perhaps the model needs tweaking.)

Plot @@ {V[t] /. First@pl, Flatten[{t, V["Domain"] /. pl}]}


We can stop the integration when V[t] == 0 with WhenEvent. Note here if you comment out Method -> {"DiscontinuityProcessing" -> False}, we get the 1/0 errors, if you use the OP's Cd. It definitely seems like bug in NDSolve.

pl = NDSolve[{eqn, bc, WhenEvent[V[t] == 0, "StopIntegration"]},
{j, V}, {t, 0, 0.01}, Method -> {"DiscontinuityProcessing" -> False}]

Plot @@ {V[t] /. First@pl, Flatten[{t, V["Domain"] /. pl}]}


• Actually that chatter is what I'm most interested in. I'll be doing a further analysis on it to see if I can obtain bifunctional or chaotic behavior by varying A and ω. I'll look this over carefully and see if I can streamline the output. The starting value of V[0] was purely a guess at what the amplitude would be, and figured it wouldn't matter after a short amount of time. Commented Jun 8, 2015 at 0:02
• @LeviJohannesHall Interesting. If you are concerned about discriminating between what behavior is due to the model and what is due to numerical truncation & round-off error, using different settings for WorkingPrecision sometimes can help sort that out. (In this example, the chatter seems well-behaved, however.) Commented Jun 8, 2015 at 4:21
• Thank you! This will definitely be of use. Commented Jun 8, 2015 at 18:25

It is a physical phenomenon occurring in a very short time( milliseconds).

A = 1;
ϕ = 0.0257275275;
ω = 151.6*10^3;
Cj = 17*10^(-9);
Cs = 35*10^(-9);
Rt = 223.9;
I0 = 5*10^-6;
L = 15*10^-3;
Cd[t_] := Piecewise[{{Cj/Sqrt[1 - A*Cos[ω*t]/ϕ],
t < 0}, {Cs*Exp[A*Cos[ω*t]/ϕ], t >= 0}}];

eqn = {V'[t] == (j[t] - I0*(Exp[A*Cos[ω*t]/ω] - 1))/
Cd[t], j'[t] == (A*Cos[ω*t] - Rt*j[t] - V[t])/L};
bc = {V[0] == A, j[0] == 0};
pl = NDSolve[{eqn, bc}, {V, j}, {t, 0, 0.0001},
MaxStepFraction -> 0.0001] // Quiet;


Plots:

Plot[V[t] /. pl, {t, 0, 0.000011}, PlotRange -> Full,
AxesLabel -> {"t[ \[Mu]s ]", "V"},
Ticks -> {{#, 10^6 #} & /@ FindDivisions[{0., 10^-5}, 10],Automatic}]
Plot[V[t] /. pl, {t, 0.000011, 0.0000132},
AxesLabel -> {"t[ \[Mu]s ]", "V"},
Ticks -> {{#, 10^6 #} & /@ FindDivisions[{0., 10^-4}, 100] // N,
Automatic}]
Plot[j[t] /. pl, {t, 0, 0.0000149}, PlotRange -> Full,
AxesLabel -> {"t[ \[Mu]s ]", "j[ m A ]"},
Ticks -> {{#, 10^6 #} & /@ FindDivisions[{0., 10^-5}, 5] //
N, {#, 10^5 #} & /@ FindDivisions[{-10^-4, 10^-4}, 5] // N}]