1
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ – Mr Puh May 4 '20 at 21:02
  • $\begingroup$ 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..... $\endgroup$ – user62716 May 4 '20 at 21:06
  • $\begingroup$ Thank you Nasser, take your time please $\endgroup$ – user62716 May 4 '20 at 22:51
  • 2
    $\begingroup$ Please put more effort in understanding answers you've obtained, there's no essential difference between your recent FDM questions. $\endgroup$ – xzczd May 5 '20 at 2:38
  • $\begingroup$ Thank you xzczd, I have got the perfect answer from Dr.Nasser. $\endgroup$ – user62716 May 5 '20 at 11:54
4
$\begingroup$

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

Mathematica graphics


enter image description here


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
 ]
$\endgroup$
6
  • $\begingroup$ 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 $\endgroup$ – user62716 May 5 '20 at 11:40
  • $\begingroup$ @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. $\endgroup$ – Nasser May 5 '20 at 16:21
  • $\begingroup$ 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. $\endgroup$ – Nasser May 5 '20 at 16:26
  • $\begingroup$ 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 $\endgroup$ – user62716 May 5 '20 at 19:21
  • $\begingroup$ @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. $\endgroup$ – Nasser May 5 '20 at 21:25

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.