5
$\begingroup$

I would like to solve two different partial differential equations each one defined in a different region and in different coordinates. However the equations are coupled through a boundary condition at a common boundary. How can I do this in Mathematica?

To simplify the problem and make the question more specific. If, for example, I consider time-dependent diffusion problem on the space interval [0, 2] (assuming in the following diffusion coefficient(s) to be equal to unity):

eqns={D[c1[t,x],t]==D[c1[t,x],x,x]};
bcs={c1[t,0]==1,c1[t,2]==Exp[-1000*t]};
ics={c1[0,x]==1};

Tmax=0.5;
sol=NDSolveValue[{eqns,bcs,ics},c1,{t,0,Tmax},{x,0,2}]
Plot3D[sol[t,x],{t,0,Tmax},{x,0,2}, PlotRange->All]

This works fine. If now I reformulate the same problem for two diffusion equations defined at the regions [0,1] and [1,2], respectively, and connect them through the equality of concentrations and fluxes at the point x=1, I would assume this should result in something like this:

eqns={D[c1[t,x],t]==D[c1[t,x],x,x],D[c2[t,y],t]==D[c2[t,y],y,y]};
bcs={c1[t,0]==1,c1[t,1]==c2[t,1],Derivative[0, 1][c1][t, 1] == Derivative[0, 1][c2][t, 1],c2[t,2]==Exp[-1000*t]};
ics={c1[0,x]==1,c2[0,y]==1};

Tmax=0.5;
sol=NDSolveValue[{eqns,bcs,ics},{c1,c2},{t,0,Tmax},{x,0,1},{y,1,2}]

The latter, however return the error "The length of the derivative operator Derivative[1,0] in c1(1,0)[t,x] is not the same as the number of arguments." assuming that both c1 and c2 should be the functions of t, x and y and not c1[t,x], c2[t,y].

Could you please help to solve the problem?

$\endgroup$

1 Answer 1

8
$\begingroup$

If I understand your question correctly, you could solve one PDE over the entire region and have different PDE coefficients active in different parts of the region, like so:

sol = NDSolveValue[{Inactive[
       Div][-If[x <= 1, {{10}}, {{1}}].Inactive[Grad][
         u[x], {x}], {x}] == 1, DirichletCondition[u[x] == 0, True]}, 
   u, {x, 0, 2}];
Plot[sol[x], {x, 0, 2}]

enter image description here

Update

As long as the point where you want an internal DirichletCondition is part of the mesh you can have those too:

sol = NDSolveValue[{Inactive[
       Div][-If[x <= 1, {{10}}, {{1}}].Inactive[Grad][
         u[x], {x}], {x}] == 1, 
    DirichletCondition[u[x] == 0, x == 0 || x == 2], 
    DirichletCondition[u[x] == 0.5, x == 1]}, u, {x, 0, 2}, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"IncludePoints" -> {{1.}}}}];
Plot[sol[x], {x, 0, 2}]

Note that everything remains the same, I just included a point at coordinate 1. and set a different DirichletCondition there.

enter image description here

$\endgroup$
6
  • $\begingroup$ Thanks for help! I believe this will work in my case. Since I need to solve different PDEs in two contiguous regions, as I stated in the question, I can combine both equations into a single general PDE and zero some terms in one region and the rest of the terms in the other, like you proposed. Thank you once again. $\endgroup$
    – Alexander
    May 22, 2016 at 10:38
  • $\begingroup$ The drawback of this approach, however, is that it's impossible to set explicitly a boundary condition at the "sewing" point. $\endgroup$
    – Alexander
    May 22, 2016 at 13:25
  • 1
    $\begingroup$ @Alexander, there is no such drawback. $\endgroup$
    – user21
    May 22, 2016 at 20:43
  • $\begingroup$ Agree! ) Thanks a lot! $\endgroup$
    – Alexander
    May 23, 2016 at 20:27
  • $\begingroup$ Just out of curiosity, what if I want a 3rd type boundary condition using NeumannValue[] at x==1? $\endgroup$
    – 407PZ
    Sep 29, 2017 at 9:15

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.