3

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[1][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[0] == 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?

xzczd
  • 65,995
  • 9
  • 163
  • 468
Lufu
  • 51
  • 2
  • Can't you just integrate your n function by hand and explicitly add it to the equation? – gpap Apr 15 '14 at 09:55
  • The actual problem is a little bit more complicated, so there is no formal solution for n[t] – Lufu Apr 15 '14 at 09:57
  • Can you post a minimal example capturing the intricacies of the original problem then? – gpap Apr 15 '14 at 09: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.8u[t, x] - (5.n[t]Derivative[0, 1][u][t, x])/ (5. + n[t]) + (50.n[t]^2Derivative[0, 2][u][t,x])/(5. + n[t])^2, Derivative[1][n][t] == 100. - 0.8n[t] - (0.5Integrate[u[t,x], {x, 0, 2Pi}]n[t])/(5. + n[t]), u[0,x] == 200., n[0] == 50., u[t, 0] == u[t, 2Pi]}, {u, n}, {t, 0., 2.}, {x, 0., 6.28}] This should be correct now. – Lufu Apr 15 '14 at 11:39
  • AFAIK NDSolve can't automatically handle this, you can couple PDEs with PDEs and ODEs with ODEs but not PDEs with ODEs. But as NDSolve for all versions up to 9 uses the method by lines only for PDEs you can just as well do the method by lines by hand: generate a set of coupled ODEs by discretizaing in space and then couple that system to the single ODE. See tutorial/NDSolvePDE in the online documentation for more details, an example and even some internal helper methods which help with the discretization... – Albert Retey Apr 15 '14 at 12:12
  • Is this related? http://mathematica.stackexchange.com/q/22189/78 – P. Fonseca Apr 15 '14 at 19:19

1 Answers1

6

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[NIntegrate`GaussRuleData[points, MachinePrecision]];    
midgrid = Rescale[nodes, {0, 1}, domain];

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

grid = Flatten[{domain[[1]], 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

Mathematica graphics

Plot[nsol[t], {t, 0, tend}, PlotRange -> All]

Mathematica graphics

xzczd
  • 65,995
  • 9
  • 163
  • 468