Re: How to specify boundary conditions on all 4 sides of a plate for a steady state heat equation (PDE) using NDSolve? (Laplace equation)
- To: mathgroup at smc.vnet.net
- Subject: [mg59780] Re: How to specify boundary conditions on all 4 sides of a plate for a steady state heat equation (PDE) using NDSolve? (Laplace equation)
- From: Peter Pein <petsie at dordos.net>
- Date: Sun, 21 Aug 2005 03:51:27 -0400 (EDT)
- References: <ddpt58$orc$1@smc.vnet.net> <ddscf6$ai1$1@smc.vnet.net> <de6m0a$ch8$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Nasser Abbasi schrieb: > "Jens-Peer Kuska" <kuska at informatik.uni-leipzig.de> wrote in message > news:ddscf6$ai1$1 at smc.vnet.net... > >>Hi, >> >>NDSolve[] is for intial-boundary value problems >>n+1 and not >>for pure boundary value problems. >>You can use the tim depend equation and integrate >>it until the solution >>does not change any more. >> >>Regards >> Jens > > > Well, I think I'll go back to using direct numerical/mesh based > methods. > > Those will work for any type of linear PDE (in general), and even if > one must write more code to do that, one can get more control on what > is going on, and one does not have the limitations imposed by NDSolve. > > For fun, here is the solution to 1-D diffusion PDE using FTCS scheme > (forward time centered space) using Mathematica code. I just hacked > this quickly, and would like to go over it to again to make it more of > a 'functional' style, but I find it hard to avoid using the For loop > sometimes, but I do use the Table command a lot, so I think the code > is functional :) > > http://12000.org/my_notes/mma_matlab_control/e65/HTML/e65.htm > > compare the above to the solution given by using NDSolve shown below > (much less code ofcourse) > > http://12000.org/my_notes/mma_matlab_control/e57/HTML/e57.htm > > Next, I'll solve the 2D steady state heat equation (Laplace PDE) using > such direct method, may be using a different scheme. > > btw, I found Mathematica Matrix operations and general performance > doing this to be really fast, I was surprised, I thought it will be a > little slower. > > Nasser > > Hi Nasser, if you want to avoid For[] loops, try NestList[] (needs a lot of memory in the line ttplot=...): Off[General::spell1]; makeD[n_] := Module[{idn = Flatten[IdentityMatrix[n]]}, idn = {1, -2, 1} . (RotateLeft[idn, #1] & ) /@ {-1, 0, 1}; idn[[{1, 2}]] = idn[[{-2, -1}]] = {0, 0}; Partition[idn, n]] nPoints = 31; len = 4; \[Kappa] = 1/2; h = len/(nPoints - 1); \[Tau] = 1/10^4; k = (\[Tau]*\[Kappa])/h^2; plotStep = (nSteps = 10^4)/(nplots = 50); tt = N[Sin[h*Pi*(Range[nPoints] - 1)], 16]; bb = N[Inverse[IdentityMatrix[nPoints] - k*makeD[nPoints]], 16]; ttplot = NestList[bb.#&, tt, nSteps][[1 + Range[0, nSteps, plotStep]]]; ListPlot3D[ttplot, MeshRange -> {{0, len}, {0, \[Tau]*nSteps}}, AxesLabel -> {"x", "time/s", "T[x,t]"}, PlotRange -> All, ImageSize -> 500]; -- Peter Pein Berlin http://people.freenet.de/Peter_Berlin/