Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

What I'm trying to achieve is model of the heat flow, in this case for the simplest 1D case,its relatively easy to do for the steady state case, but when I try to do it with NDsolve so I get the distribution of heat over time,I fail to come up with a good way to connect the two differential equaions, I recently modeled a gravity field where the same differential equation is applied many times to different bodies, that I got to work, but I dont understand how to link "differeet" ODE/PDE's

This is my code, so far, and I intend to display that solution with Manipulate and Plot:

:ps I added the [] in u[1] and u[2] last minute to see if it changes anything before they were just called u1[] and u2[] respectively, it might be nonsensical but my original question remains

thanks in advance!

heateqs = {D[u[1][x, t], t] == D[u[1][x, t], {x, 2}], 
   D[u[2][x, t], t] == 0.5 D[u[2][x, t], {x, 2}]};

bcs = {u[1][0, t] == 300, u[1][5, t] == u[2][5, t], 
   Derivative[1, 0][u[2]][10, t] == 0};

ics = {u[1][x, 0] == u[2][x, 0] == 0};

sol = NDSolve[{heateqs, bcs, ics}, u, {x, 0, 10}, {t, 0, 10}]

NDSolve::bcedge: Boundary condition u[1][5,t]==u[2][5,t] is not specified on a single edge of the boundary of the computational domain. >>
share|improve this question
why you have two equations for 1D problem? isn't this is a 1D heat equation? it should be only one equation which works for both time and space. –  Algohi Jun 25 '14 at 1:08
well if I assume that the first half of my 1dimensional cable is 1/2 L long and made out of say copper and the second half is made from Iron, they will conduct heat at different rates, following different PDE's but the heat flux entering one would be the same as exiting the other,and since the heat flux at the border is dependtant on temperatures it cant be solved in independant eq's either it seems –  catadoxas Jun 25 '14 at 1:19

2 Answers 2

I think you don't need to have two equation to describe behavior of domain follows same PED. you may need to add UnitStep to the thermal diffusivity factor.

check if this work for you (domain of length=1):

sol = u[x, t] /. NDSolve[{
(1 - 0.5 UnitStep[x - 0.5]) D[u[x, t], t] == D[u[x, t], x, x],
u[x, 0] == 0,
u[0, t] == 300,
    (D[u[x, t], x] /. x -> 1) == 0
   u, {t, 0, 10}, {x, 0, 1}]

With[{sol = sol},
 Manipulate[Plot[sol, {x, 0, 1}, PlotRange -> {0, 350}], {t, 0, 10}]]

you may notice that there is boundary and initial condition inconsistency. you have u=0 when t=0 & any x and at the same time you have u=300 at x=0 and any t.

share|improve this answer
I don't think this describes the physics expressed by the first two conditions properly. –  Jens Jun 25 '14 at 2:54
why not? the thermal conductivity is constant over each half of the domain. but I general your solution is much better than this solution. –  Algohi Jun 25 '14 at 2:56
Well, I'm assuming that the heat source is at x=0 and therefore we need to have u=300 there at all times... but probably it's best to let the OP decide what conditions they want to impose. –  Jens Jun 25 '14 at 2:58
actually in the real case the value of u is not just jump from 0 to 300 one shot. instead, it goes gradually(some cases rapidly). that is why we do have inconstancy in B.C and I.C. In general you give the solution for catadoxas and I think he will like your solution. –  Algohi Jun 25 '14 at 3:03
OK, I'll upvote your solution because really the basic idea we use is identical, only the physical interpretation differs! –  Jens Jun 25 '14 at 3:08

You don't need two different functions because the position-dependent material parameter can be incorporated into a function that I'll call d[x] and that enters in a single heat conduction equation as follows:

d[x_] := (1 + 4 UnitStep[5 - x])/5.

heateq = d[x] D[u[x, t], t] == D[u[x, t], {x, 2}];

Plot[d[x], {x, 0, 10}, PlotRange -> {0, 1}, Exclusions -> None]


Here, I exaggerated the jump from large to small thermal conductivity so that the interface becomes more easily visible in the plot of the solution below.

To solve the differential equation, we need to add one more condition to the two main constraints that you require (starting height at x=0 is 300, and spatial derivative vanishes at x=10). I added the functional form at t=0 in the form of a narrow Gaussian:

sol = First[NDSolve[{heateq,
     u[x, 0] == 300 Exp[-x^2 10],
     u[0, t] == 300,
     Derivative[1, 0][u][10, t] == 0.}, 
    u[x, t], {x, 0, 10}, {t, 0, 10.}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "DifferenceOrder" -> "Pseudospectral"}}]

frames = Table[
   Plot[Evaluate[u[x, t] /. sol], {x, 0, 10}, 
    PlotRange -> {0, 300}], {t, .001, 10, .2}];

ListAnimate[Show[#, Graphics[Line[{{5, 0}, {5, 300}}]]] & /@ frames]


The vertical line indicates the interface where the jump in d[x] occurs. In the Method option of NDSolve, I chose the method of lines, following the suggestion in this answer.

share|improve this answer
I get similar idea like you and I was working on it. when I post me solution I found your solution. maybe I need to delete me solution which is definitely not as good as your solution . –  Algohi Jun 25 '14 at 2:53
Thank you for the answer. This is a very interesting problem. However, I think that flux continuity is not satisfied: -(1.0 Derivative[1, 0][u][5.0, t]) should be equal to -(1/5 Derivative[1, 0][u][5.0, t]). What do you think? –  Laurent Jan 13 at 12:44
@Laurent In order to get a jump in the derivative, one would have to specify an interface condition that is clearly not present in the original question. With the information in the question (which is not a steady-state problem and contains no conditions on the derivative at the interface), I think we can't get what you're expecting. –  Jens Jan 13 at 16:21
Just a sidenote, I believe that OP and @Laurent and this OP did want to add a flux continuity condition, they just didn't understand what a flux continuity condition is very well. Flux continuity condition depends on heat conductivity coefficient (usually notated with $\lambda$) rather than thermal diffusivity (usually notated with $\alpha$, for OP's question, (1 + 4 UnitStep[5 - x])/5. plays the role of thermal diffusivity exactly. ) –  xzczd Apr 14 at 13:41
BTW, let alone the vagueness of this question, how to Implement heat flux continuity condition is indeed a interesting question, I do find a solution for this, but not that general and I'm not satisfied. I'm hesitating whether I should post it as an answer or start a new question, which will of course state the problem in a clearer way. –  xzczd Apr 14 at 13:46

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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