2

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.

xzczd
  • 65,995
  • 9
  • 163
  • 468
S. Maths
  • 203
  • 1
  • 9
  • 1
    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. – user21 Jan 28 '19 at 07:36
  • 1
    @S.Cho Where did this problem come from? – Alex Trounev Jan 28 '19 at 13:32
  • 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. – S. Maths Jan 28 '19 at 18:17
  • @user21 I edited my first post. – S. Maths Jan 28 '19 at 22:13
  • 1
    Related: https://mathematica.stackexchange.com/q/144881/1871 – xzczd Feb 01 '19 at 15:57
  • Thank you ! This helped me to understand many new things to simplify the problem. It still a smal problem with NDSolve as follows https://i.stack.imgur.com/ndRR0.jpg – S. Maths Feb 01 '19 at 22:05
  • I think this clarify the first one https://i.stack.imgur.com/OjSbC.jpg – S. Maths Feb 01 '19 at 23:11

2 Answers2

3

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

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106
3

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • Thank you Sir. Very interesting !. It seems to work in a good way. I need some time to better understand the pdetoode function. Is It works when the domain is a Disk for example ? Thanks again. – S. Maths Jan 29 '19 at 11:49
  • 1
    @S.Cho It can handle any regular domain (rectangular, disk, etc.). – xzczd Jan 29 '19 at 12:03
  • @xzczd Is it possible to expand the use of the function pdetoode to solve 3D and 4D? – Alex Trounev Jan 30 '19 at 12:38
  • 1
    @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. – xzczd Jan 30 '19 at 12:51
  • @xzczd I tried with centered initial condition u[0,x]==Exp[-(x-0.5)^2] but the plotting not working. What's the problem ? – S. Maths Jan 30 '19 at 13:27
  • @xzczd Thank you. I reproduced the example. I will try to apply pdetoode for other problems. – Alex Trounev Jan 30 '19 at 13:34
  • 1
    @S.Cho Try adding PlotRange->All to Plot3D. – xzczd Jan 30 '19 at 13:58
  • @xzczd I tried to Solve the 2D case by pdetoode (the same equation in Disk[] using polar cordinate, since the derivative on the boundary are not easy with Cartesian cordinate) but It doesn't work. Should I change the pdetoode to fix this ? – S. Maths Jan 30 '19 at 19:49
  • In fact I tried this one

    `domain = Disk[]; 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, r, z], 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, r, z], {t, 0, tend}, {r, 0, 1}, {z, 0, 2 Pi}, PlotRange -> All] & @@ domain`

     – S. Maths Jan 30 '19 at 19:49
  • with this `With[{u = u[t, r, z]}, {eq, ic, bc} = {D[u, t] == D[u, r, r] + (1/r) D[u, r] + (1/r^2) D[u, z, z], u == 1 /. t -> 0, {(D[u, t] == D[u, r, r] + (1 + 1/r) D[u, r] + (1/r^2) D[u, z, z]) /. r -> 1}}];

    test = u -> (1 &);

    {eq, ic, bc} /. test`

     – S. Maths Jan 30 '19 at 19:52
  • Sorry in this case the domain with r and z is a Rectangle. but still this problem. https://i.stack.imgur.com/Cjvm0.jpg – S. Maths Jan 30 '19 at 21:03
  • 1
    @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, #&. – xzczd Jan 31 '19 at 12:38
  • Thank you for suggestions. I know the method of lines in 1-D, here is 2-D problem, and in fact I'm new to using mathematica. I tried to Solve the problem basing on the previous example but I don't get where is the problem. For example what are the points deleted here xdomain = ydomain = {-1, 1}; del = Delete[#, {{1}, {2}, {-2}, {-1}}] &; – S. Maths Jan 31 '19 at 21:13
  • @S.Cho Notice pdetoode generates ODEs for every grid points when discretizing PDE, while we also have equations generated from b.c., so we need to remove some of the equations to make room for the b.c.. If you still don't get it, try a coarse grid and observe the output carefully. – xzczd Feb 01 '19 at 13:10
  • The problem still in how to choose these equations to be removed. First I thought that it is related to the boundary of the domain, but It seems not. – S. Maths Feb 02 '19 at 17:24
  • 1
    @S.Cho You may want to read the comments under this post. – xzczd Feb 03 '19 at 02:26
  • Still misundesrtanding the problem, may be because I'm new to Mathematica and also my english is not good enough to understand better. I think a simple example in 2-D should help. I posted my problem here https://mathematica.stackexchange.com/questions/190754/ndsolve-error-in-2-d-heat-equation. Thank you for every things. – S. Maths Feb 03 '19 at 12:39
  • 1
    @S.Cho One thing I forgot to say, "I know the method of lines in 1-D, here is 2-D problem"，then I'm afraid you don't know method of lines. Method of lines for 1D and 2D have no essential difference. Once again, pdetoode is just a auxiliary tool for the implementation of method of lines, you need to know the basic of this method. To test if you've understood the basic: do you know how to implement the method without pdetoode? – xzczd Feb 04 '19 at 02:35
  • What I know is : In the method of lines we descretize the equation and BCs w.r.t space variables and solving the resulting ODEs, I implemented it in 1-D. It is the same for 2-D case, but still some troubles in removing some equations to make room for bc equations as you said. If you recommend any reference to better understanding MOL It would be nice. Thank you again. – S. Maths Feb 04 '19 at 08:53
  • @S.Cho Then how will you implement method of lines without pdetoode? Notice I need to remove equations because pdetoode has been designed to generate equations for every grid point (it uses one-sided formula at and near the boundary, if you're not familiar with one-sided formula, start from page 6 of this book. This book is also a good reference for FDM, which is closely related to method of lines). – xzczd Feb 04 '19 at 10:07