NDSolve two PDE plus one ODE, together

I have a couple of PDE (time and space/length) that are correctly managed by NDSolve, but to which I’m adding an ODE on which the acc “accumulates”, through time, a “value” coming from one of the boundaries of the PDE (the b*h[t, le]*u[t, le]).

On my following example, the accumulation is not linked, and so acc could be calculated afterwards. But the purpose is to use the acc[t] value as a definition of one of the boundary of h. Once the errors are gone, I'll define something like h[t, le]==acc/100 instead of the current fixed value boundary.

I’ve tried a lot of assorted stuff, but NDSolve keeps giving errors on the amount of parameters of the ODE, etc, as shown below.

I’m now thinking on managing the “accumulation” with the EvaluationMonitoring, but it doesn’t sound right...

a = 10;
le = 62;
g = 9.8;
tf = 100;
b = 3.4;
r = 5.86/le;
NDSolve[{
b*u[t, x]*Derivative[0, 1][h][t, x] +
b*h[t, x]*Derivative[0, 1][u][t, x] +
b*Derivative[1, 0][h][t, x] == r,

(r*u[t, x])/(b*g*h[t, x]) + (Abs[u[t, x]]*(b + 2*h[t, x])^(4/3)*
u[t, x])/(a^2*(b*h[t, x])^(4/3)) +
Derivative[0, 1][h][t, x] + (u[t, x]*Derivative[0, 1][u][t, x])/
g + Derivative[1, 0][u][t, x]/g == 0,

Derivative[acc][t] == b*h[t, le]*u[t, le],

h[0, x] == 1,
u[0, x] == 0,
h[t, le] == 1,
u[t, 0] == 0,
acc == 0
},
{u, h, acc},
{t, 0, tf}, {x, 0, le}
];

(*NDSolve::ndode: Input is not an ordinary differential equation. >>*)

EDIT

I came across a problem with the naming of the variable acc. But since it is probably related with a different question, I wrote my problem here

• As there is no direct dependence of PDE to the function acc[t], can not it be solved later as separate ODE. Mar 27 '13 at 13:52
• @PlatoManiac as stated on the problem description, the simplified example doesn't have a dependence. But once the errors are gone, I'll add it to the h[t, le] boundary (actually, something like h[t, le]==acc/100). I just thought one problem at a time. (acc is a volume, 100 is an area, and h is a height...) Mar 27 '13 at 13:53
• @RolfMertig but it handles well the two coupled PDE (the result looks correct when cross-checking with other sources). It's only when I add the ODE that it gets confused. Mar 27 '13 at 15:45
• why not do Derivative[1, 0][acc][t, x] == b*h[t, le]*u[t, le] and adjust the boundary conditions Mar 27 '13 at 16:14
• @RolfMertig I get NDSolve::delpde: Delay partial differential equations are not currently supported by NDSolve. >>. Since there is no "delayed" problem, I imagine that there's a different way do define it. I've done: Derivative[1, 0][acc][t, x] == b*h[t, le]*u[t, le], acc[0, x] == 0 Mar 27 '13 at 16:20

not a answer - just a elements of information - reliability not garanted

The first problem one encounters with your code is : "NDSolve::derlen: The length of the derivative operator Derivative1 in (acc^[Prime])[t] is not the same as the number of arguments". That's not the main problem.

To get rid of it, I have created a function acc[]` of 2 variables : t and x, and stated the problem as 3 PDEs.

Derivative[acc][t] == b*h[t, le]*u[t, le] becomes :

Derivative[1, 0][acc][t, x] == b*h[t, x]*u[t, x]

Initial value acc[t]==0 becomes : acc[0,x]==0

The solution of interest is acc[t,le].

Here is the code :

NDSolve[{
b*u[t, x]*Derivative[0, 1][h][t, x] +
b*h[t, x]*Derivative[0, 1][u][t, x] +
b*Derivative[1, 0][h][t, x] == r,

(r*u[t, x])/(b*g*h[t, x]) + (Abs[u[t, x]]*(b + 2*h[t, x])^(4/3)*
u[t, x])/(a^2*(b*h[t, x])^(4/3)) +
Derivative[0, 1][h][t, x] + (u[t, x]*Derivative[0, 1][u][t, x])/
g + Derivative[1, 0][u][t, x]/g == 0,

Derivative[1, 0][acc][t, x] == b*h[t, x]*u[t, x],

h[0, x] == 1,
u[0, x] == 0,
h[t, le] == 1,
u[t, 0] == 0,
acc[0, x] == 0}, {u, h, acc}, {t, 0, tf}, {x, 0, le}]

Mathematica accepts this and gives this solution : I obtain the same with your code.

• It already helps a little. I'm checking if results are correct. But the definition of acc for a range of x is a little artificial, considering that acc is a single point in space... But thank you anyway. Mar 27 '13 at 16:31
• After checking, it returns something very wrong. Probably beacuase we are relating u and h in all their x dimension with acc, when acc is actually only influencing the le point. Mar 27 '13 at 16:43
• @P.Fonseca I have checked too and it seems not to be wrong. See the graphic in the answer. Is that what you expect ? Agree that the method is far from ideal. I investigate because I'm interested. Mar 27 '13 at 20:01
• It toke me some time to figure out the problem. Take a look at my following answer: mathematica.stackexchange.com/questions/15672/… Mar 28 '13 at 9:16
• P. Fonseca Thanks. Otherwise, my Mathematica version in 8.0.4. Mar 28 '13 at 13:41