Skip to main content

numerical integration - Incorrect results of diffusion equation with Neumann boundary conditions



I want to resolve a PDE model, which is 1D heat diffusion equation with Neumann boundary conditions. The key problem is that I have some trouble in solving the equation numerically. Consider the following code:


h = 6000;

a = 200;
Dif = 3.67*10^-14*10^18;
Ni = 1;
deq = D[u[t, x], t] == Dif*D[u[t, x], {x, 2}]
ic = u[0, x] == If[0 <= x <= a , Ni, 0]
bc = {Derivative[0, 1][u][t, 0] == 0, Derivative[0, 1][u][t, h] == 0}
sol = NDSolve[{deq, ic, bc}, u, {t, 0, 60}, {x, 0, h}]
Plot3D[Evaluate[u[t, x] /. sol], {t, 0, 60}, {x, 0, h}, PlotStyle -> Automatic]

enter image description here



I got a result, but a error was occurred.


NDSolve::ibcinc:

I know that this error suggests conflicts between initial condition and boundary conditions, although I have no idea where conflict come from.


In addition, as you can see, the value of x=0 is gradually increased with time in spite of Neumann conditions.


Any suggestions how to fix it?



Answer



The comment of @xzczd is very pertinent, but there are a lot of things to say about this subject. Among theses things :





  • In your example NDSolve automatically chooses the "TensorProductGrid" method (as opposed to "FiniteElement"). This information is sometimes hard to find. I get it from experience (Edit here is a question that asks how to know which method NDSolve has automatically chosen).




  • This choice leads to the problem mentionned by @xzczd. This problem is complicated to analyse and it is not clearly documented. I'm speaking of this documentation




  • A more friendly approach is to use the Finite Element Method. With this method, the syntax for the Neumann boundary condition is not Derivative[0, 1][u][t, 0] == 0 but a syntax that use NeumannValue. The use of NeumannValue is a little bit disturbing at the beginning, but in your case it's very simple because the boundary condition equivalent to Derivative[0, 1][u][t, 0] == 0 is the default choice of NDSolve with the finite element method.




So, to get the solution, just remove the boundary conditions and impose the finite elemnt method :



h = 6000;
a = 200;
Dif = 3.67*10^-14*10^18;
Ni = 1;
deq = D[u[t, x], t] == Dif*D[u[t, x], {x, 2}]
ic = u[0, x] == If[0 <= x <= a , Ni, 0]
sol = NDSolve[{deq, ic}, u, {t, 0, 60}, {x, 0, h},Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]
Plot3D[Evaluate[u[t, x] /. sol], {t, 0, 60}, {x, 0, h}, PlotStyle -> Automatic]


enter image description here


Comments

Popular posts from this blog

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1....