I am trying to solve this 1-D heat equation in the interval [-L1,L2] with the follow conditions:thermal diffusivity=alpha1 between [-L1,0);thermal diffusivity=alpha2 between (0,L2].Boundary conditions: T=T0+T1*Sin(omega*t) at x=-L1 and x=L2; dT/dx=flux at x=0. What's the correct way to set this up? Thanks.

  • 2
    $\begingroup$ Just play around a bit with NDSolve. $\endgroup$ Commented Jul 24, 2018 at 7:10
  • 3
    $\begingroup$ between [L1,0) you mean [-L1,0] and then you say Initial condition: T=T0+T1*Sin(omega*t) how could initial conditions have t in it? This makes it then T=T0, since sin(0)=0 because t=0 at initial conditions. As for solving, you can try to make the thermal diffusivity as Piecewise and see if NDSolve accepts it. So use BC for the whole bar as is, and plugin for alpha the Piecwise condition. $\endgroup$
    – Nasser
    Commented Jul 24, 2018 at 7:18
  • 1
    $\begingroup$ The FEM backend of NDSolve can definately treat discontinuous diffusion constants. @Nasser $\endgroup$ Commented Jul 24, 2018 at 9:15
  • $\begingroup$ I made some corrections to the question. $\endgroup$ Commented Jul 24, 2018 at 18:02
  • $\begingroup$ You've removed "Initial condition" from the question, but this doesn't clarify anything, what's the "T=T0+T1*Sin(omega*t)"? Also, what have you tried? $\endgroup$
    – xzczd
    Commented Jul 25, 2018 at 4:16

1 Answer 1


The 1D heat equation in x and t

pde = alpha[x]*D[u[x, t], x, x] - D[u[x, t], t] == 0

Its easiest for me to set up alpha as a function of UnitStep

alpha[x_] = alpha1 (UnitStep[x + L1] - UnitStep[x]) + alpha2 (UnitStep[x] - UnitStep[x - L2])

The boundary conditions.

bc1 = u[-L1, t] == T0 + T1 Sin[omega t]

bc2 = u[L2, t] == T0 + T1 Sin[omega t]

You did not have an initial condition and MMA didn't complain when I didn't add one, but I added one anyway.

ic = u[x, 0] == T0

Plug in some numbers

alpha1 = 1;
alpha2 = 2;
T0 = 1;
T1 = 2;
omega = 2;
L1 = 1;
L2 = 2;

NDSolve[{pde, bc1, bc2, ic}, u[x, t], {x, -L1, L2}, {t, 0, 20}] // Flatten;

u[x_, t_] = u[x, t] /. %

tp = Table[Plot[u[x, t], {x, -L1, L2}, PlotRange -> {-1, 3}], {t, 0, 20, .1}];

enter image description here


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.