# How to solve the heat flow PDE for the constant source term

I have the two dimensional heat flow equation solved for a two dimensional reagion over time. That is an instance I found in the internet. What I am looking for is a method how to determine the heat that is introduced to obtain the given temperature distribution.

In a first step I made up a code to obtain a temperature distribution with a known heat source term at given location. In a second step I take the obtained temperature distribution to obtain the source term. In this way I can test the method of obtaining the source term, before I use it in a real case.

This is the model to get the temperature distribtion for a source term:

Ω =
ImplicitRegion[! ((x - 5)^2 + (y - 5)^2 <= 3^2), {{x, 0, 5}, {y, 0,
10}}];
op = \!$$\*SubscriptBox[\(∂$$, $$t$$]$$u[t, x, y]$$\) - \!$$\*SubsuperscriptBox[\(∇$$, $${x, y}$$, $$2$$]$$u[t, x, y]$$\);
Subscript[Γ, D] =
DirichletCondition[u[t, x, y] == 0, x == 0 && 8 <= y <= 10];
Subscript[Γ, N] =
NeumannValue[-1*u[t, x, y], (-5 + x)^2 + (-5 + y)^2 == 9] +
NeumannValue[10(*source term to be determined later again*), 2 <= x <= 5 && y == 0];
uifHeat =
NDSolveValue[{op == Subscript[Γ, N],
Subscript[Γ, D], u[0, x, y] == 0},
u, {t, 0, 100}, {x, y} ∈ Ω];
Manipulate[
ContourPlot[
uifHeat[t, x, y], {x, y} ∈ uifHeat["ElementMesh"],
ColorFunction -> own, Contours -> 20, AspectRatio -> 2,
PlotRange -> All], {{t, 10}, 0, 100, 2}]
f[x_] = uifHeat[100, x, 10];(*This will be used as input for second step*)


Sorry it is here a bit hard to read but copying it into Mathematica formats it much more readable.

In the second step I added the temperature distribution f[x] as a DirichletCondition and replace the 10 in the NeumannValue by a constant P. The way tried it gives no solution but errors:

Argument {{{{P}},{{P}},{{P....}}}} at position 1 should be a rank 3 tensor of machine-size real numbers

The code is:

op = \!$$\*SubscriptBox[\(∂$$, $$t$$]$$u[t, x, y]$$\) - \!$$\*SubsuperscriptBox[\(∇$$, $${x, y}$$, $$2$$]$$u[t, x, y]$$\);
Subscript[Γ,
D] = {DirichletCondition[u[t, x, y] == 0, x == 0 && 8 <= y <= 10],
DirichletCondition[u[100, x, y] == f[x],
0. <= x <= 5. && y == 10.]};
Subscript[Γ, N] =
NeumannValue[-1*u[t, x, y], (-5 + x)^2 + (-5 + y)^2 == 9] +
NeumannValue[P (*To be determined by NDSovleValue*), 2 <= x <= 5 && y == 0];
{uifHeat, h} =
NDSolveValue[{op == Subscript[Γ, N],
Subscript[Γ, D], u[0, x, y] == 0}, {u, T}, {t, 0,
100}, {x, y} ∈ Ω];


Can I determine P without trial and error? How do I do it?

****Edit 01.09.15: It seems not clear what the task is. I am looking for a way to determine the heating power necessary to for a temperature profile that was/will be measured (indirect power determination by simulation). The first block of code is only to get a termperature profile f[x] under a known condition of heating (NeumannValue[10,...]) and cooling. The second code block should contain a procedure to calculate/simulate the heating power when f[x], as an additional boundary condition, is given and the source term NeumannValue[10,...] only contains a free paremeter P (NeumannValue[P,..]. P should evaluate to app. 10 within the procedure I'm looking for, as this is the value for the Neumann condition when f[x] was calculated with. I figured out a solution meanwhile but it works extremly slow. This is the code

ClearAll[lsg]
lsg[P_Real] :=
NDSolveValue[{op ==
NeumannValue[-1*u[t, x, y], (-5 + x)^2 + (-5 + y)^2 == 9] +
NeumannValue[P, 2 <= x <= 5 && y == 0],
Subscript[\[CapitalGamma], D], u[0, x, y] == 0},
u, {t, 0, 100}, {x, y} \[Element] \[CapitalOmega]]
AbsoluteTiming[NMinimize[
{dd = lsg[P];
Total[Table[dd[100, x, 10 ] - f[x], {x, 0, 5, 0.1}]^2], P > 9.}
, P, MaxIterations -> 2, Method -> "NelderMead"]]


The output is:

{517.944, {7.66488*10^-17, {P -> 10.}}} I get the error "Failed to converege..." but still the solution is right.

• If you're solving a PDE numerically, then T needs to have a specific numerical value. There's really no way around it. – Michael Seifert Aug 27 '15 at 13:30
• So it can not be solved or does it just need a starting value? – Eisbär Aug 27 '15 at 14:21
• I'm honestly not 100% sure what you're after, but if I understand you correctly: no, NDSolve can't find the conditions on one part of the boundary that lead to a particular solution at the end of the time interval. If I've misunderstood your question, though, you might want to edit it to make it clearer. – Michael Seifert Aug 27 '15 at 19:09
• One should be able to use ParametricNDSolve to create ParametricFunction representing the solution as a function of T. This could be used to solve for T using an appropriate solver. – Michael E2 Aug 28 '15 at 10:31
• @Mr.Wizard Classes just started, so my free time is somewhat random. When this first appeared, I didn't have time, simply left a hint, and forgot all about it. Like Michael Seifert, I'm not sure of what is sought. Perhaps it's in the code somewhere how to determine whether one has obtained "the given temperature distribution," but I don't see it. If I knew that, I probably could answer. As for closing, I would say it's not a "simple mistake" but perhaps "unclear" (needs more info). It might turn out to be a duplicate if the question were clearer (haven't looked, though). – Michael E2 Aug 31 '15 at 22:58