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

sol@Vs = rebuild[uval[[2]], 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 '18 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 '18 at 21:09
• I realized how to plot the boundary variables. – kosa Sep 27 '18 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[1][[3 ;; 4]]` at $x=0$ and $x=1$, so we need to make room for them. – xzczd Sep 28 '18 at 6:01
• Yes, makes sense. I will ptoofunc, to understand the spatial discretization. – kosa Sep 28 '18 at 6:44