# NDSolve: Couple PDE and ODE involving integral

I am trying to numerically solve an equation with NDSolve, where there is a ODE coupled to a PDE, like the following:

    NDSolve[{
Derivative[1, 0][u][t,x] ==
-0.8*u[t, x] - (5.*n[t]*Derivative[0, 1][u][t, x])/(5. + n[t]) + (50.*n[t]^2*Derivative[0, 2][u][t, x])/(5. + n[t])^2,
Derivative[n][t] ==
100. - 0.8*n[t] - (0.5*Integrate[u[t, x], {x, 0, 2*Pi}]*n[t])/(5. + n[t]),
u[0, x] == 200.,
n == 50.,
u[t, 0] == u[t, 2*Pi]},
{u, n},
{t, 0., 2.},
{x, 0., 6.28}]


Unfortunately Mathematica tells me "NDSolve::ndode: Input is not an ordinary differential equation".

When i change n[t] to n[t,x] it calculates something, but n[t,x] doesn't stay uniform over time, which it should because n wasn't a function of x in the first place.

Would anyone know a way around this?

• Can't you just integrate your n function by hand and explicitly add it to the equation? – gpap Apr 15 '14 at 9:55
• The actual problem is a little bit more complicated, so there is no formal solution for n[t] – Lufu Apr 15 '14 at 9:57
• Can you post a minimal example capturing the intricacies of the original problem then? – gpap Apr 15 '14 at 9:58
• Sorry I can't make sense of this as there are a bunch of brackets missing, all the variables have been switched to [Theta], [Omega] in place of \[Theta], \[Omega] and a constant isn't defined. You could edit the question with the full problem - that would help people get interested in this. – gpap Apr 15 '14 at 10:17
• NDSolve[{Derivative[1, 0][u][t,x] == -0.8*u[t, x] - (5.*n[t]*Derivative[0, 1][u][t, x])/ (5. + n[t]) + (50.*n[t]^2*Derivative[0, 2][u][t,x])/(5. + n[t])^2, Derivative[n][t] == 100. - 0.8*n[t] - (0.5*Integrate[u[t,x], {x, 0, 2*Pi}]*n[t])/(5. + n[t]), u[0,x] == 200., n == 50., u[t, 0] == u[t, 2*Pi]}, {u, n}, {t, 0., 2.}, {x, 0., 6.28}] This should be correct now. – Lufu Apr 15 '14 at 11:39

## 1 Answer

As mentioned in the comment above, NDSolve doesn't know how to disretize your system, so we do it ourselves. I'll use pdetoode for the disretization of the derivative term and Gaussian quadrature formula for the discretization of integral:

With[{u = u[t, x], n = n[t]},
eq@1 = D[u, t] == -0.8 u - (5. n D[u, x])/(5. + n) + (50. n^2 D[u, x, x])/(5. + n)^2;
eq@2 = D[n, t] == 100. - 0.8 n - (0.5 (Integrate[u, {x, 0, 2 Pi}]) n)/(5. + n);
{ic@1, ic@2} = {u == 200, n == 5} /. t -> 0];

points = 25; difforder = 4; domain = {0, 2 Pi};

{nodes, weights} = Most[NIntegrateGaussRuleData[points, MachinePrecision]];
midgrid = Rescale[nodes, {0, 1}, domain];

intrule = HoldPattern@
Integrate[__] -> -Subtract @@ domain weights.Map[u[#][t] &, midgrid];

grid = Flatten[{domain[], midgrid, domain[[-1]]}];

(*Definition of pdetoode isn't included in this post,
please find it in the link above.*)
ptoofunc = pdetoode[u[t, x], t, grid, difforder, True];
ode@1 = ptoofunc@eq@1;
ode@2 = eq@2 /. intrule;
{odeic@1, odeic@2} = {ptoofunc@ic@1, ic@2};

tend = 2;
{usollst, nsol} =
NDSolveValue[{ode /@ {1, 2}, odeic /@ {1, 2}}, {u /@ grid, n}, {t, 0, tend}];
usol = rebuild[usollst, grid];


Notice the periodic b.c. is set with the 5th argument of pdetoode.

Plot3D[usol[t, x], {t, 0, tend}, {x, ##}] & @@ domain Plot[nsol[t], {t, 0, tend}, PlotRange -> All]
` 