2
$\begingroup$

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.

$\endgroup$
  • 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 '19 at 7:36
  • 1
    $\begingroup$ @S.Cho Where did this problem come from? $\endgroup$ – Alex Trounev Jan 28 '19 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 '19 at 18:17
  • $\begingroup$ @user21 I edited my first post. $\endgroup$ – S. Maths Jan 28 '19 at 22:13
  • 1
    $\begingroup$ Related: mathematica.stackexchange.com/q/144881/1871 $\endgroup$ – xzczd Feb 1 '19 at 15:57
3
$\begingroup$

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"}]}

fig1

|improve this answer|||||
$\endgroup$
  • 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 '19 at 9:53
  • 1
    $\begingroup$ @S.Cho The method is called explicit Euler. $\endgroup$ – Alex Trounev Jan 29 '19 at 11:56
  • 2
    $\begingroup$ @xzczd Thank you! I updated the message. $\endgroup$ – Alex Trounev Jan 29 '19 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 '19 at 14:05
  • 1
    $\begingroup$ @S.Cho Time is lost on building ListPointPlot3D[] $\endgroup$ – Alex Trounev Jan 29 '19 at 15:15
3
$\begingroup$

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

$$u(t,x)=1$$

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

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ @S.Cho It can handle any regular domain (rectangular, disk, etc.). $\endgroup$ – xzczd Jan 29 '19 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 '19 at 12:51
  • 1
    $\begingroup$ @S.Cho Try adding PlotRange->All to Plot3D. $\endgroup$ – xzczd Jan 30 '19 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 '19 at 12:38
  • 1
    $\begingroup$ @S.Cho You may want to read the comments under this post. $\endgroup$ – xzczd Feb 3 '19 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.