Solving a second order coupled PDE system

Question

This is similar to How to solve a certain coupled first order PDE system but I seem to be getting errors which is most likely due to my misunderstanding on how the code is actually working. (I'm really bad at coding)

I have taken 2 different approaches to the problem, one is using the method from the link above, the other is using code I wrote. The problem with the linked code, is that I do not know how to add my boundary and initial conditions.

The system I would like to solve is as follows, where $U=U(x,t);V=V(x,t)$

$a \frac{\partial U}{\partial t} = b \frac{\partial^2 U}{\partial x^2} +cV \exp(-\frac{1}{U})$

$d \frac{\partial V}{\partial t} + e\frac{\partial V}{\partial x} = f \frac{\partial^2 V}{\partial x^2} - g V \exp(-\frac{1}{U})$

where $a-g$ and $U_0,V_0$ are constants to be decided later.

The i.cs and b.cs for the system are:

$x=0 : \frac{\partial V(0,t)}{\partial x}=0 , \frac{\partial U(0,t)}{\partial x}=0 $

$x=1 : \frac{\partial V(1,t)}{\partial x}= V-V_0 , \frac{\partial U(1,t)}{\partial x}=U-U_0 $

$t=0 : V(x,0)=V_0,U(x,0)=U_0$

Firstly, My code is as follows:

ClearAll["Global'*"]
System = {a D[U[x, t], t] == b D[U[x, t], x, x] + c V[x, t]*Exp[-1/U[x, t]],
  d D[V[x, t], t] + e D[V[x, t], x] == 
   f D[V[x, t], x, x] - g V[x, t]*Exp[-1/U[x, t]], D[V[0, t], x] == 0, 
  D[V[1, t], x] == V[x,t]-v, D[U[0, t], x] == 0, D[U[1, t], x] == U[x,t]-u, 
  V[x, 0] == v, U[x, 0] == u}
NDSolve[System, {U[x, t], V[x, t]}, {x, 0, 5}, {t, 0, 5}];

Plotting wont be much of an issue as instead of a manipulate plot, those constants will just be given values in the actual equations later.

The second code is as from the link here on stack exchange

r0 = First@
   Solve[{z[x, t] == u[x, t] - v[x, t], 
     y[x, t] == u[x, t] + v[x, t] }, {u[x, t], v[x, t]}];
eq1 = (Unevaluated[
      a D[u[x, t], t] - b D[u[x, t], x, x] - c v[x, t] + 
       Exp[1/u[x, t]]] /. %) // Simplify;
eq2 = (Unevaluated[
      d D[v[x, t], t] + e D[v[x, t], x] - f D[v[x, t], x, x] - 
       g v[x, t] Exp[1/ u[x, t]]] /. %%) // Simplify;
Simplify[eq1 - eq2];
r1 = Simplify[#] & /@ DSolve[% == 0, z[x, t], {x, t}][[1, 1]]
Simplify[eq1 + eq2 ];
r2 = Simplify[#] & /@ DSolve[% == 0, y[x, t], {x, t}][[1, 1]] /. 
  C[1] -> C[2]
r0 /. {r1, r2}

Not only does the above code confuse me, but I have no idea how I will plot the functions nor do I know how to incorporate my boundary conditions. But the main issue with this code is I don't understand the mathematica commands and method used, so ideally my code (if possible) would be better to edit.

Any help would be greatly appreciated

Many thanks

Ken

Answer 1

Issue 1: You get the error

NDSolve : Equation or list of equations expected instead of True in the first argument

because, e.g., D[V[1, t], x] returns 0 and System contains equations of the form D[V[1, t], x] == 0 which evaluates to True. See, V[1, t] does really not depend on x. What you mean is D[V[x, t], x] == 0 /. x -> 0 where you may read /. x -> 0 as "at $x = 0$".

The cleaned up System looks like this:

System = {
  a D[U[x, t], t] == b D[U[x, t], x, x] + c V[x, t] Exp[-1/U[x, t]], 
  d D[V[x, t], t] + e D[V[x, t], x] == f D[V[x, t], x, x] - g V[x, t] Exp[-1/U[x, t]],
  D[V[x, t], x] == 0 /. x -> 0,
  D[V[x, t], x] == V[x, t] - v /. x -> 1,
  D[U[x, t], x] == 0 /. x -> 0,
  D[U[x, t], x] == U[x, t] - u /. x -> 1,
  V[x, 0] == v,
  U[x, 0] == u
  }

Issue 2: It should be something like

  NDSolve[System, {U[x, t], V[x, t]}, {x, 0, 1}, {t, 0, 5}];

since otherwise, the spacial boundary is not the boundary.

Issue 3: Using NDSolve with symbols appearing in the equations is futile. But after setting them to numeric values, it works just fine:

ClearAll[U, V]
a = 1; b = 1; c = 1; d = 1; e = 1; f = 1; g = 1;
v = 6; u = 4;
System = {
   a D[U[x, t], t] == b D[U[x, t], x, x] + c V[x, t] Exp[-1/U[x, t]], 
   d D[V[x, t], t] + e D[V[x, t], x] == f D[V[x, t], x, x] - g V[x, t] Exp[-1/U[x, t]],
   D[V[x, t], x] == 0 /. x -> 0,
   D[V[x, t], x] == V[x, t] - v /. x -> 1,
   D[U[x, t], x] == 0 /. x -> 0,
   D[U[x, t], x] == U[x, t] - u /. x -> 1,
   V[x, 0] == v,
   U[x, 0] == u
   };
{U, V} = NDSolveValue[System, {U, V}, {x, 0, 1}, {t, 0, 5}]

Plot3D[{U[x, t], V[x, t]}, {x, 0, 1}, {t, 0, 5}]