# Both DSolve and NDSolve fail for a simple system of two PDEs

Using Mathematica 10.4.1, I'm trying to solve some simple (I think!) partial differential equations with a form more or less like this:

$$\qquad \frac{\partial v}{\partial t} = b v\;, \quad \frac{\partial v}{\partial p} = k v$$

with the initial condition $$v(t_0,p_0)=v_0$$. My actual problems of interest are a bit more complex, but if I can't get this simple case to work then anything more complicated is hopeless.

I'm not a specialist in partial differential equations, but I'd think this would be a relatively straightforward problem: by eye, it certainly seems like $$v = v_0 e^{b(t-t_0)+k(p-p_0)}$$ is the solution. But when I try to use DSolve on this problem:

DSolve[
{D[v[t, p], t] == b v[t, p], D[v[t, p], p] == k v[t, p], v[t0, p0] == v0},
v[t, p], {t, p}]


Mathematica just spits my input back at me. I'd really like to think that a system this simple would be within Mathematica's capacity to handle! Am I doing something wrong?

Understanding that closed-form PDE solutions might just not be viable, I figured I'd just see about getting a numerical solution, so I replaced the constants with arbitrary numbers and tried what I thought would be sure to give a brute-force numerical integration of the system:

NDSolve[
{D[v[t, p], t] == 2 v[t, p], D[v[t, p], p] == 3 v[t, p], v[300, 1] == 5},
v, {t, 100, 500}, {p, 0.1, 5}]


But instead, I get an error:

NDSolve::overdet: "There are fewer dependent variables, {v[t,p]}, than equations, so the system is overdetermined."

I can't for the life of me understand how this is overdetermined. I've got a function of two variables: shouldn't specifying both partials and an initial condition be exactly what you need for a PDE? What am I missing here?

• It is overdetermined, but that does not mean there are no special solutions. Try solving each PDE separately and see if you can reconcile them. Apr 29, 2019 at 15:47
• If this is a system of PDE, then initial and boundary conditions must be specified. Condition v[300, 1] == 5 does not correspond to initial or boundary conditions. Apr 29, 2019 at 17:45
• @MichaelE2: Given that the mixed partials commute, I had thought there was a theorem somewhere stating that a solution should exist. Again: not my specialty! Maybe there's some limitation on the form of the functions that specify the partials that I'm not aware of. (But your overall point is well taken: treating the two variables separately seems like a sensible approach, and in this case it does reproduce the solution that I guessed by eye. I'm not sure how to apply that when seeking a numerical solution, though.) Apr 29, 2019 at 19:07
• @AlexTrounev: I'm not sure what sort of initial/boundary condition you have in mind as more appropriate here. As I said originally, this particular system has a straightforward solution, and that solution includes only a single undetermined constant. Specifying the value of the function at one point is sufficient to uniquely determine that constant. I feel like the same would be true for a broad class of systems of this type: roughly speaking, I imagine numerically integrating to find v(t,p0) for all t, and then integrating along p to fill the region. Is there a flaw in that intuition? Apr 29, 2019 at 19:15
• @Steuard PDE is defined over the region, in this case on a rectangle {t, 100, 500}, {p, 0.1, 5} . Depending on the type of equations, a certain number of conditions are set at the boundary of the region. But not at one point inside the region. Apr 29, 2019 at 22:15

This is a system of hyperbolic equations. Differentiating the first equation with respect to p, the second with respect to t, then, since $$\frac {\partial ^2v}{\partial p\partial t}=\frac {\partial ^2v}{\partial t\partial p}$$, we have $$b\frac {\partial v}{\partial p}=k\frac {\partial v}{\partial t}$$ Or simply take v from the first equation and substitute it into the second.For this equation it is necessary to formulate a problem in a rectangular region. The original system of equations can be used to find the boundary conditions, for example, when p = 0.1, we get from the first equation v[t,.1]==v0*Exp[b*t], and at t=100 we have v[100,p]==v1*Exp[k*p]. But numerically we will not be able to solve this problem, since there are given large numbers like Exp. A solution using ODE is discussed in the topic Numerical methods to solve a continuity equation . A simple example of an analytical solution.

b = 2; k = 3; v0 = 1; v1 = v0; sol =
DSolve[{b*D[v[t, p], p] == k*D[v[t, p], t], v[t, 0] == v0*Exp[b*t],
v[0, p] == v1*Exp[k*p]}, v, {t, 0, 5}, {p, 0, 5}]

(*Out[]= {{v -> Function[{t, p}, E^(3 p + 2 t)]}}*)

Plot3D[v[t, p] /. sol, {t, 0, 5}, {p, 0, 5}, Mesh -> None,
PlotRange -> All] A simple example of a numerical solution.

b = 2; k = 3; v0 = 1; v1 = v0; sol =
NDSolve[{b*D[v[t, p], p] == k*D[v[t, p], t], v[t, 0] == v0*Exp[b*t],
v[0, p] == v1*Exp[k*p]}, v, {t, 0, 5}, {p, 0, 5}]

Plot3D[v[t, p] /. sol, {t, 0, 5}, {p, 0, 5}, Mesh -> None,
PlotRange -> All]

• Thank you! For the benefit of others with similar questions, here are the two big messages that I've taken from this answer (and from your comments earlier). First, while my single data point was sufficient to make the problem well-posed, Mathematica's DE solvers need you to translate that into actual boundary data. And second, equality of mixed partials means that my two DEs aren't as independent as they look, so combining them as you have here isn't "throwing away information" (or, not much). And the info that was lost wound up encoded in the extended boundary data (coming full circle). Apr 30, 2019 at 12:47
• @Steuard On the toy model, everything is clear, but it would be better if you describe your real problem. Apr 30, 2019 at 14:32
• Now that I know how to make Mathematica do what I need for the toy model, I figure I have a good start on knowing how to make headway myself on my real questions. (My philosophy is usually to ask for help on minimal relevant test cases, and to spare others the effort of dealing with calculations that I'm able to make headway on myself. Also, I hope that the lessons of the toy model will generalize more easily for future readers, who might get lost in the weeds of a more involved system. But I don't know how well that matches the usual culture here.) Apr 30, 2019 at 18:27

This can be solved one equation at a time.

eqns = {D[v[t, p], t] == b v[t, p],
D[v[t, p], p] == k v[t, p], v[t0, p0] == v0};


First equation

sol1 = DSolve[First[eqns], v, {t, p}]
(* {{v -> Function[{t, p}, E^(b t) C[p]]}} *)

eqns2 = eqns /. First[sol1]
(* {True, E^(b t) Derivative[C][p] == E^(b t) k C[p],
E^(b t0) C[p0] == v0} *)


Second equation

sol2 = DSolve[eqns2[], C, p]
(* {{C -> Function[{p}, E^(k p) C]}} *)

eqn3 = eqns2 /. First[sol2]
(* {True, True, E^(k p0 + b t0) C == v0} *)


Third equation

sol3 = Solve[eqn3[], C]
(* {{C -> E^(-k p0 - b t0) v0}}


Then we can put it back to get the full result.

v[t, p] /. First[sol1] /. First[sol2] /. First[sol3] *)
(* E^(k p - k p0 + b t - b t0) v0 *)


I imagine that there are some sets of equations where this approach would give an answer that is incorrect or incomplete.