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

d[x]

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]

diffuse

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 at 2:53

Your Answer

 
discard

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.