Finite difference method for 1D heat equation

I have solve the following 1D heat equation: ut=uxx, t>0,0<=x<=5, with ic=u(x,0)=x^2, and bcs:u(0,t)=2t;u(5,t)=2t+25:

ClearAll["Global*"]
heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == (x^2); bc = {u[0, t] == 2 t, u[5, t] == 2 t + 25};
sol = DSolve[{heqn, ic, bc }, u[x, t], {x, t}]


I got the exact solution u=2t+x^2.

Now I implement the explicit finite difference method: ut=ui,j+1-ui,j/delta t; uxx=ui-1,j-2ui,j +ui+1,j/delta x, then I got:

ui,j+1=r[ui-1,j-2ui,j+ ui+1,j]+ui,j; where r=delta t/delta x; then:

ui,j+1=rui-1,j+(1-2r)ui,j+rui+1,j;.................(1)

Now using Eq.(1), I want to construct A, b and getting u at x=1,2,3,4; t=0.25,0.5,0.75,1, I do not know how to get the required values numerically?

• could you rephrase the question you have and format the numerical part maybe? As it is right now it is hard to understand your question. May 4, 2020 at 21:02
• Dear Mr Puh, the question is simply, apply the finite difference method for 1D heat equation, the formulations used for ut, uxx are given, we need to find u at some points at given time values..... May 4, 2020 at 21:06
• Thank you Nasser, take your time please May 4, 2020 at 22:51
• Please put more effort in understanding answers you've obtained, there's no essential difference between your recent FDM questions. May 5, 2020 at 2:38
• Thank you xzczd, I have got the perfect answer from Dr.Nasser. May 5, 2020 at 11:54

This is basic explicit method finite difference. An implicit method will be better and left as an exercise.

This is described in https://en.wikipedia.org/wiki/Finite_difference_method#Explicit_method

makeA[n_] := Module[{A, i, j},
A = Table[0, {i, n}, {j, n}];
Do[
Do[
A[[i, j]] = If[i == j, -2, If[i == j + 1 || i == j - 1, 1, 0]],
{j, 1, n}
],
{i, 1, n}
];
A[[1, 1]] = 1;
A[[1, 2]] = 0;
A[[-1, -1]] = 1;
A[[-1, -2]] = 0;
A];

makeInitialU[nPoints_, L_, h_, ic_, leftBC_, rightBC_] :=
Module[{u, j, t = 0},
u = Table[0, {j, nPoints}];
Do[
u[[j]] = If[j == 1, leftBC[0, 0],
If[j == nPoints, rightBC[L, 0],
ic[(j - 1)*h, 0]]
],
{j, 1, nPoints}
];
u
];

updateU [currentU_, currentTime_, nPoints_, L_, h_, initialC_,
leftBC_, rightBC_, delT_, diffusion_, A_] := Module[{u},
u = ((delT/h^2) * diffusion*A . currentU) + currentU;
u[[1]] = leftBC[0, currentTime]; (*set to BC condition*)
u[[-1]] = rightBC[L, currentTime];(*set to BC condition*)
u
];

ic[x_, t_] := x^2;
leftBC[x_, t_] := 2 t;
rightBC[x_, t_] := 2 t + 25;

Manipulate[
Module[{nextU , g, currentTime = 0, j, currentU, L, h, A,
exactSolution, pExact, pFDM, k, x, t, tmp, last},

exactSolution = 2 maxtime + x^2;
L = 5;
h = L/(nPoints - 1);
currentU = makeInitialU[nPoints, L, h, ic, leftBC, rightBC];
A = makeA[nPoints];
(*iteration loop to update FDM in time*)
Do[
currentTime = currentTime + delT;
last = currentU;
currentU = updateU[currentU, currentTime, nPoints, L, h, ic, leftBC, rightBC,
delT, diffusion, A]
, {j, 0, Round[maxtime/delT]}
];

pFDM = ListLinePlot[Transpose[{Range[0, L, h], currentU}] ,
Mesh -> All, PlotStyle -> Red];

pExact = Plot[exactSolution, {x, 0, L},
AxesOrigin -> {0, 0},
PlotStyle -> Blue, ImageSize -> 400,
GridLines -> Automatic, GridLinesStyle -> LightGray,
AxesLabel -> {"x", "u(x,t)"},
PlotLabel ->
Style["Explicit method finite difference method. exact (blue) vs. FDM (red)", 12]
];
last = MatrixForm[ (NumberForm[#, {5, 2}] &) /@ last];
tmp = MatrixForm[ (NumberForm[#, {5, 2}] &) /@ currentU];

g = Grid[{{Row[{"time ", NumberForm[currentTime, {4, 2}]}], SpanFromLeft },
{Row[{" h = ", NumberForm[N@h, {5, 5}]}], SpanFromLeft },
{Row[{"CFL conditions (make sure to keep below 1/2 =  k*delt/h^2  = ",
NumberForm[diffusion*delT/h^2, {5, 5}]}], SpanFromLeft },
{Row[{"U = ", NumberForm[(delT/h^2) * diffusion, {6, 5}], " * ",
MatrixForm[A], " . ", last , " + ", last , " = ", tmp}],
SpanFromLeft},
{ Show[pExact, pFDM, PlotRange -> All], SpanFromLeft}}, Frame -> All];
g
]
,
{{nPoints, 3, "points"}, 3, 20, 1, Appearance -> "Labeled"},
{{delT, 0.01, "time space (delT)"}, 0.001, 0.1, .001, Appearance -> "Labeled"},
{{maxtime, 6, "maxtime"}, 0, 20, .1, Appearance -> "Labeled"},
{{diffusion, 1}, None},

TrackedSymbols :> {nPoints, maxtime, delT},
ContinuousAction -> False
]
`
• Great Dr.Nasser as always your work perfect and you have high level skills, please what about for Neumann and mixed boundary conditions, can you help me in this and one more please I have post 2D Laplace with mixed and no one can help me in this post....Best regards May 5, 2020 at 11:40
• @user62716 Please find more complete demonstrations on this topic I found at Wolfram demonstration site SolvingTheConvectionDiffusionEquationIn1DUsingFiniteDifferen/ and [demonstrations.wolfram.com/… these support different boundary conditions and different configurations. You can download these with the source code as well. All done using finite difference methods. May 5, 2020 at 16:21
• For the Laplace in 2D, please see this one at Wolfram demonstration site SolvingThe2DPoissonPDEByEightDifferentMethods/ it has all the options you are asking for. It also uses Finite difference method and other methods which you can choose. May 5, 2020 at 16:26
• Dear Dr.Nasser, many thanks for your help and very useful sites. Just final two things:1) can get the results not only the plots, I mean matrix, b, u in that sites?2) I will post final problem for 1D wave then will close the FDM....Best regards May 5, 2020 at 19:21
• @user62716 I think those demos at Wolfram site use Finite Difference. There are more if you search that site. I do not remember now if these ones print the A,b matrices. But you could download the code, and add a Print statement at the right place to see these matrices. For more details, I saw this page which describes Neumann Boundary conditions on 2D grid may be that will be of help. May 5, 2020 at 21:25