6
$\begingroup$

I am trying to simulate a combination of PDEs and ODEs, given below.

$$ \begin{matrix} -L\dfrac{\partial}{\partial t}I(t,z)&=&\dfrac{\partial}{\partial z}V(t,z)+RI(t,z)\\ C\dfrac{\partial}{\partial t}V(t,z)&=& -GV(t,z)-\dfrac{\partial}{\partial z}I(t,z) \end{matrix}\hspace{8mm}0<z<1 $$

$$ \begin{matrix} -L_0\dfrac{d}{d t}I(t,0)&=&R_0I(t,0)+V(t,0)-E_0(t)\\ \dfrac{d}{dt}E_0(t)&=&-K_p\left(I(t,0)-I_d\right)-K_d \dfrac{\partial}{\partial t}I(t,0) \end{matrix}~\;\;\;z=0\\ $$

$$ \begin{matrix} -L_1\dfrac{d}{d t}I(t,1)&=& R_1I(t,1)-V(t,1)+v(t)\\ C_1\dfrac{d}{dt}v(t)&=& I(t,1)-\dfrac{v(t)}{R_C} \end{matrix}~\;\;\;\;z=1 $$

All constants can be taken as 1.

PS: I have tried the following:

l = 1; l0 = 1; l1 = 1; c = 1; c1 = 1;
r = 1; g = 1; r0 = 1; r1 = 1; rc = 1;
Ki = 1;
Kd = 2;
Ib0d = 5;
simtime = 10;
var = {Is[t, x], Vs[t, x], Ib0[t], Ib1[t], Vb1[t], Eb0[t]} // Flatten
varInit = var /. {t -> 0};
pde = Inverse[
    DiagonalMatrix[{-l, c}]].{D[Vs[t, x], x] + 
     r*Is[t, x], -g*Vs[t, x] - D[Is[t, x], x]};
MatrixForm[pde];
ode = Inverse[
    DiagonalMatrix[{-l0, -l1, c1}]].{r0*Ib0[t] + Vs[t, 0] - Eb0[t], 
    r1*Ib1[t] + Vb1[t] - Vs[t, 1], Ib1[t] - Vb1[t]/rc};
MatrixForm[ode];
uode = {-Ki*(Ib0[t] - Ib0d) - Kd*(r0*Ib0[t] + Vs[t, 0] - Eb0[t])/-l0};
MatrixForm[uode];
f = {pde, ode, uode} // Flatten;
DAE = {Thread[D[var, t] == f], Ib0[t] == Is[t, 0], Ib1[t] == Is[t, 1],
    Thread[ varInit == Flatten[RandomReal[{0, 1}, {1, 6}]]]} // 
  Flatten; MatrixForm[DAE]
uval = NDSolve[DAE, {x, 0, 1}, {t, 0, simtime}]

but it is not working, giving me an underdetermined system error.

$\endgroup$
2
  • 1
    $\begingroup$ The last line of code should be uval = NDSolve[DAE, var, {x, 0, 1}, {t, 0, simtime}], but this reveals a more fundamental problem, mixing ODEs and PDEs. $\endgroup$
    – bbgodfrey
    Commented Sep 27, 2018 at 12:47
  • $\begingroup$ Yes, I will edit the mistake. But this does NOT solve the problem as @bbgodfrey pointed out. $\endgroup$
    – kosa
    Commented Sep 27, 2018 at 12:50

1 Answer 1

8
$\begingroup$

This problem can be solved with the help of pdetoode.

We first eliminate Ib0 and Ib1 from the code because DAE system is troublesome compared to ODE system for NDSolve. It's better to avoid introducing them from the very beginning to make the code cleaner, but here I'll just modify your DAE for simplicity:

newdae = DeleteCases[
  DAE /. Rule @@@ (DAE[[7 ;; 8]] /. (a_ == b_) :> (Head@a == (Function[t, #] &@b))), 
  True | Is[__Integer] == _]

Then we extract those equations involving Is and Vs:

{eq@1, eq@2} = GatherBy[newdae, FreeQ[#, Is | Vs[__]] &]

Finally, discretize those involving Is and Vs and solve the resulting system:

domain = {0, 1}; points = 25; difforder = 2;
grid = Array[# &, points, domain];
ptoofunc = pdetoode[{Is, Vs}[t, x], t, grid, difforder];
delete = #[[2 ;; -2]] &;

(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
eqlst = MapAt[delete, ptoofunc@eq@1, {1}]

uval = NDSolveValue[{eqlst, eq@2}, {Is /@ grid, Vs /@ grid, Head /@ var[[5 ;;]]}, {t, 0, 
   simtime}]

Is and Vs is discretized in $x$ direction in uval. If you prefer a continuous version, we can:

sol@Is = rebuild[uval[[1]], grid]

sol@Vs = rebuild[uval[[2]], grid]

Plot3D[sol[Vs][t, x], {t, 0, 10}, {x, 0, 1}]

Mathematica graphics

$\endgroup$
5
  • $\begingroup$ Thanks a Lot! Awesome. I still don't understand everything. I am trying to understand the code. I quick question: How do I plot the boundary variables (`Eb0[t], Vb1[t] ') $\endgroup$
    – kosa
    Commented Sep 27, 2018 at 19:33
  • $\begingroup$ in the line " eqlst = MapAt[delete, ptoofunc@eq@1, {1}]", Why are we removing the first and last lines? $\endgroup$
    – kosa
    Commented Sep 27, 2018 at 21:09
  • $\begingroup$ I realized how to plot the boundary variables. $\endgroup$
    – kosa
    Commented Sep 27, 2018 at 21:20
  • 1
    $\begingroup$ @kosa It's because ptoofunc generates discretized equation on every grid points i.e. $x=0,\frac{1}{24},\frac{2}{24},…, \frac{23}{24}, 1$, but we already have eq[1][[3 ;; 4]] at $x=0$ and $x=1$, so we need to make room for them. $\endgroup$
    – xzczd
    Commented Sep 28, 2018 at 6:01
  • $\begingroup$ Yes, makes sense. I will ptoofunc, to understand the spatial discretization. $\endgroup$
    – kosa
    Commented Sep 28, 2018 at 6:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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