3
$\begingroup$

I would like to Investigate if the edges of this 2D heat equation are insulated by the Mathematica assumed Neumann conditions.

aap = NDSolve[{\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(f[x, y, t]\)\) == 0.05 \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(f[x, y, t]\)\), 
    f[x, y, 0] == E^(-40 (x^2 + y^2))}, 
   f, {t, 0, 1}, {x, y} \[Element] Disk[]][[1]]

The problem is that Mathematica seems unable to Integrate the solution to the problem.

    Export["test.gif", 
 Table[Plot3D[aap[x, y, i], {x, y} \[Element] Disk[], 
   PlotRange -> {Automatic, Automatic, {0, 2}}], {i, .1, 5, 0.2}], 
 "AnimationRepetitions" -> \[Infinity], "ImageSize" -> 20]

enter image description here

I try

NIntegrate[f[x, y, 0] /. aap, {x, y} \[Element] Disk[]]

But my kernel just quits, same thing on Mathematica Online.

Do you also see this?/Is someone able to help?

edit:

I now did get the result one time:

In[5]:= NIntegrate[f[x, y, 0] /. aap, {x, y} \[Element] Disk[]]

During evaluation of In[5]:= Outer::heads: Heads Alternatives and List at positions 3 and 2 are expected to be the same.

During evaluation of In[5]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[5]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.07854482887754174` and 2.6180493354023447`*^-7 for the integral and error estimates.

Out[5]= 0.0785448

which is the result I need. After that I managed to crash the kernel again and I am back at the same point...

edit (2) the workaround solution as mention below works! The following now shows there is indeed no outflow of "heat".

aap = NDSolveValue[{D[f[x, y, t], t] == 
     0.05 Laplacian[f[x, y, t], {x, y}], 
    f[x, y, 0] == 40/\[Pi] E^(-40*(x^2 + y^2))}, 
   f, {t, 0, 10}, {x, y} \[Element] Disk[]];

and plot

Plot[NIntegrate[
  aap[x, y, t], {x, y} \[Element] aap["ElementMesh"]], {t, 0, 10}, 
 PlotRange -> {0, 2}]

to get:

enter image description here

$\endgroup$
5
$\begingroup$

That's a bug which I filed. As a workaround you can use:

aap = NDSolveValue[{D[f[x, y, t], t] == 
     0.05 Laplacian[f[x, y, t], {x, y}], 
    f[x, y, 0] == E^(-40*(x^2 + y^2))}, 
   f, {t, 0, 1}, {x, y} \[Element] Disk[]];
NIntegrate[aap[x, y, 0], {x, y} \[Element] aap["ElementMesh"]]
0.07854480753557473`
$\endgroup$

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.