# Simulating a combination of PDEs and ODEs

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

$$\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.

• 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. Sep 27, 2018 at 12:47
• Yes, I will edit the mistake. But this does NOT solve the problem as @bbgodfrey pointed out.
– kosa
Sep 27, 2018 at 12:50

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,
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[], grid]

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

Plot3D[sol[Vs][t, x], {t, 0, 10}, {x, 0, 1}] • 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] ')
– kosa
Sep 27, 2018 at 19:33
• in the line " eqlst = MapAt[delete, ptoofunc@eq@1, {1}]", Why are we removing the first and last lines?
– kosa
Sep 27, 2018 at 21:09
• I realized how to plot the boundary variables.
– kosa
Sep 27, 2018 at 21:20
• @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[[3 ;; 4]]` at $x=0$ and $x=1$, so we need to make room for them. Sep 28, 2018 at 6:01
• Yes, makes sense. I will ptoofunc, to understand the spatial discretization.
– kosa
Sep 28, 2018 at 6:44