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.

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

2 Answers 2

Reset to default
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$
13
  • 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 ♦
    Commented Jan 29, 2019 at 9:53
  • 1
    $\begingroup$ @S.Cho The method is called explicit Euler. $\endgroup$
    – Alex Trounev
    Commented Jan 29, 2019 at 11:56
  • 2
    $\begingroup$ @xzczd Thank you! I updated the message. $\endgroup$
    – Alex Trounev
    Commented Jan 29, 2019 at 12:07
  • 1
    $\begingroup$ @S.Cho It is necessary to compare the two methods on non-trivial initial data. $\endgroup$
    – Alex Trounev
    Commented Jan 29, 2019 at 14:05
  • 1
    $\begingroup$ @S.Cho Time is lost on building ListPointPlot3D[] $\endgroup$
    – Alex Trounev
    Commented Jan 29, 2019 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$
20
  • 1
    $\begingroup$ @S.Cho It can handle any regular domain (rectangular, disk, etc.). $\endgroup$
    – xzczd ♦
    Commented Jan 29, 2019 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 ♦
    Commented Jan 30, 2019 at 12:51
  • 1
    $\begingroup$ @S.Cho Try adding PlotRange->All to Plot3D. $\endgroup$
    – xzczd ♦
    Commented Jan 30, 2019 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 ♦
    Commented Jan 31, 2019 at 12:38
  • 1
    $\begingroup$ @S.Cho You may want to read the comments under this post. $\endgroup$
    – xzczd ♦
    Commented Feb 3, 2019 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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