Can we solve the following PDE by Mathematica,

enter image description here

where $\Omega$ is a bounded domain of $\mathbb{R}^n$, $\Gamma =\partial \Omega$ is the boundary of $\Omega$, $\partial_\nu$ is the normal derivative, and $\nu$ is the outer unit vector.

First I tried to solve the equation by Matlab, but the tools (pdepe, pde toolbox.. etc) not include dynamic boundary conditions $(2)$.

I tried the NDSolve with mathematica, but this not work for me:

NDSolve[{D[u[t, x], t] == D[u[t, x], {x, 2}], u[0, x] == 1, 
Derivative[1, 0][u][t, 0] == Derivative[0, 1][u][t, 0],
Derivative[1, 0][u][t, 1] == -Derivative[0, 1][u][t, 1]}, u, {t, 0, 10}, {x, 0, 1}]

It generate the following error:

NDSolve::bdord: Boundary condition -(u^(0,1))[t,0]+(u^(1,0))[t,0] should have derivatives of order lower than the differential order of the partial differential equation.

  • 1
    $\begingroup$ Have you tried anything? If yes, what issues did you encounter? It would make it much easier for people to answer your question if you provided the mathematica code you have so that we do not have to type it. $\endgroup$ – user21 Jan 28 at 7:36
  • 1
    $\begingroup$ @S.Cho Where did this problem come from? $\endgroup$ – Alex Trounev Jan 28 at 13:32
  • $\begingroup$ This is a heat equation with dynamic boundary conditions which involves time derivative on the boundary, Yes I tried many tools in Matlab as in Mathematica (NDSolve...) but all the known methods don't work for this type of boundary conditions. $\endgroup$ – S. Maths Jan 28 at 18:17
  • $\begingroup$ @user21 I edited my first post. $\endgroup$ – S. Maths Jan 28 at 22:13
  • 1
    $\begingroup$ Related: mathematica.stackexchange.com/q/144881/1871 $\endgroup$ – xzczd Feb 1 at 15:57

We can build a converging process.

U[0][x_] := 1; t0 = 1/50;
Do[U[t] = 
   NDSolveValue[{(u[x] - U[t - t0][x])/t0 == 
      D[u[x], x, x], (u[0] - U[t - t0][0])/t0 == 
      u'[0], (u[1] - U[t - t0][1])/t0 ==- u'[1]}, u, {x, 0, 1}];, {t, 
  t0, 10, t0}]
lst = Table[{t, x, U[t][x]}, {t, 0, 10, t0}, {x, 0, 1, .02}];
{ListPointPlot3D[lst, AxesLabel -> {"t", "x", "u"}, PlotRange -> All],
  Plot[U[10][x], {x, 0, 1}, AxesLabel -> {"x", "u"}]}


  • 1
    $\begingroup$ 1.Your 2nd b.c. is different from OP's. 2. For the b.c. you've chosen, your solution is not reliable, it's easy to verify the correct solution is $u(t,x)=1$. $\endgroup$ – xzczd Jan 29 at 9:53
  • 1
    $\begingroup$ @S.Cho The method is called explicit Euler. $\endgroup$ – Alex Trounev Jan 29 at 11:56
  • 2
    $\begingroup$ @xzczd Thank you! I updated the message. $\endgroup$ – Alex Trounev Jan 29 at 12:07
  • 1
    $\begingroup$ @S.Cho It is necessary to compare the two methods on non-trivial initial data. $\endgroup$ – Alex Trounev Jan 29 at 14:05
  • 1
    $\begingroup$ @S.Cho Time is lost on building ListPointPlot3D[] $\endgroup$ – Alex Trounev Jan 29 at 15:15

First of all it's not hard to notice the solution for your toy example is


With[{u = u[t, x]},
  {eq, ic, bc} =
   {D[u, t] == D[u, x, x],
    u == 1 /. t -> 0,
    {D[u, t] == D[u, x] /. x -> 0, D[u, t] == -D[u, x] /. x -> 1}}];

test = u -> (1 &);

{eq, ic, bc} /. test
(* {True, True, {True, True}} *)

Anyway, it's enough for illustrating a solution so I won't change the example. Here I'd like to show a solution by discretizing the PDE in $x$ direction ourselves. I'll use pdetoode to facilitate the generation of ODE system:

domain = {0, 1};
difforder = 4;
points = 25;
grid = Array[# &, points, domain];
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[u[t, x], t, grid, difforder];
removeredundant = #[[2 ;; -2]] &;
ode = removeredundant@ptoofunc@eq;
odebc = ptoofunc@bc;
odeic = ptoofunc@ic;
tend = 10;
sollst = NDSolveValue[{ode, odeic, odebc}, u /@ grid, {t, 0, tend}];
sol = rebuild[sollst, grid]

Plot3D[sol[t, x], {t, 0, tend}, {x, ##}] & @@ domain

Mathematica graphics

  • 1
    $\begingroup$ @S.Cho It can handle any regular domain (rectangular, disk, etc.). $\endgroup$ – xzczd Jan 29 at 12:03
  • 1
    $\begingroup$ @AlexTrounev pdetoode /pdetoae can handle high dimensional PDEs (here is a post using pdetoode to solve a 2+1D problem), but the troublesome parts are: 1. Problem that is spatially 3D or higher is not that easy to solve for any numeric solver in most cases, AFAIK; 2. Additionally, since pdetoode/pdetoae generates symbolic ODE/AE for every grid points, the memory requirement of pdetoode/pdetoae is rather demanding when the problem is spatially 3D or higher, which is usually beyond the reach of PC. $\endgroup$ – xzczd Jan 30 at 12:51
  • 1
    $\begingroup$ @S.Cho Try adding PlotRange->All to Plot3D. $\endgroup$ – xzczd Jan 30 at 13:58
  • 1
    $\begingroup$ @S.Cho You're using pdetoode blindly, notice it's just a auxiliary tool for the implementation of method of lines, so you need to know the basic of method of lines. Also, though you don't need to read through the source code of pdetoode, you must understand what those codes outside of pdetoode are doing, including but not limited to Array, #&. $\endgroup$ – xzczd Jan 31 at 12:38
  • 1
    $\begingroup$ @S.Cho You may want to read the comments under this post. $\endgroup$ – xzczd Feb 3 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.