How to specify boundary conditions on all 4 sides of a plate for a steady state heat equation (PDE) using NDSolve? (Laplace equation)
Nasser Abbasi
hi;
just for fun, I am trying to solve a steady state heat equation i.e.
laplace equation, for a rectangular plate.
So, I have 4 boundary conditions, one for each side of the plate.
But when I do that, NDSolve says that it is designed to solve initial
conditions problems only? is this really the case? May be I am not
defining the B.C. correctly for Mathematica?
The code is below, also I've posted it on my web page with the full
error message.
http://12000.org/my_notes/mma_matlab_control/e61/HTML/e61.htm
I find the error strange, saying that NDSolve can only solve IC PDE,
because I solved 1-D heat equation using IC and BC earlier with no
problem, see this
http://12000.org/my_notes/mma_matlab_control/e57/HTML/e57.htm
So, I have a feeling that NDSolve can do this, I must be just doing
something not right.
Remove["Global`*"];
h = 30; w = 10; temp = 100;
eq = D[T[x, y], x, x] + D[T[x, y], y, y] == 0;
bc = {T[0, y] == 0, T[w, y] == 0, T[x, 0] == temp,T[x, h] == 0};
sol = NDSolve[{eq, bc}, T[x, y], {x, 0, w}, {y, 0, h}]
"Boundary values may only be specified for one independent
variable. Initial values may only be specified at one value of the
other independent variable."
Nasser
James Gilmore
This issue is dealt with in the help browser _explicitly_: NDSolve->Further
Examples->P.D.Es: Limitations
This is a classic mathematical physics BVP. You should approach this problem
in Mathematica, as you would by hand: use separation of variables, and then
a fourier expansion to satisfy the boundary conditions. There are many books
that explain how to do this. I look forward to tackling a similar problem by
hand when I sit my electrodynamics qualifying exam in 3 weeks.
James Gilmore
Graduate Student
Department of Physics
Yale University
New Haven, CT 06520 USA
Email: james....@yale.edu
URL: <http://pantheon.yale.edu/~jbg39/>
"Nasser Abbasi" <n...@12000.org> wrote in message
news:ddpt58$orc$1...@smc.vnet.net...
Jens-Peer Kuska
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
"Nasser Abbasi" <n...@12000.org> schrieb im
Newsbeitrag news:ddpt58$orc$1...@smc.vnet.net...
Nasser Abbasi
news:dds9co$91c$1...@smc.vnet.net...
> This is a classic mathematical physics BVP. You should approach this
> problem
> in Mathematica, as you would by hand: use separation of variables,
> and then
> a fourier expansion to satisfy the boundary conditions.
Not if you want to use a numerical solvers such as NDSolve. That is
the whole idea of using NDSolve.
I know how to solve these by hand, and also by direct numerical
approach, I've solved my of these before and more advanced ones when I
took some courses at the Math dept at UC Berkeley one year ago, I was
just playing around to see if NDSolve can solve BVP and get the same
plots I got when I solved this problem by hand using sepration of
variables.
> There are many books
> that explain how to do this.
Yes, and my home library contains many fine books on PDE's. I like the
Satnley Farlow book, and Mary Boas has excellent chapter on the
subject, but a bit short on detailes. Also Richard Haberman applied
PDE's is nice, and if you want to see a nice new book with the cover
showing solutions of PDE's plots which I am sure was made using
Mathematica check Charles MacCluer's BVP and fourier expansions Dover
book, but it does not contain any Mathematica code. Another book which
uses a CAS system to solve PDE's is by David Betounes called PDE's for
computational science, lots of examples and plots.
Is there a short list somewhere which makes it clear what kind/class
of PDE's Mathematica can solve and not solve directly using NDSolve or
even DSolve? And why is it that NDSolve can solve an initial value PDE
and not BVP? I wonder if NDSolve will be able to solve a BVP PDE in
next version?
Ferdinand Cap
F. Cap, Mathematical Methodes in Physics and Engineering with
Mathematica, CRC PRess, 2003, ISBN 1-58488-402-9
On Wed, 17 Aug 2005, Nasser Abbasi wrote:
> "James Gilmore" <james....@yale.edu> wrote in message
> news:dds9co$91c$1...@smc.vnet.net...
>
>> This is a classic mathematical physics BVP. You should approach this
>> problem
>> in Mathematica, as you would by hand: use separation of variables,
>> and then
>> a fourier expansion to satisfy the boundary conditions.
>
> Not if you want to use a numerical solvers such as NDSolve. That is
> the whole idea of using NDSolve.
>
> I know how to solve these by hand, and also by direct numerical
> approach, I've solved my of these before and more advanced ones when I
> took some courses at the Math dept at UC Berkeley one year ago, I was
> just playing around to see if NDSolve can solve BVP and get the same
> plots I got when I solved this problem by hand using sepration of
> variables.
>
>> There are many books
>> that explain how to do this.
>
> Yes, and my home library contains many fine books on PDE's. I like the
> Satnley Farlow book, and Mary Boas has excellent chapter on the
> subject, but a bit short on detailes. Also Richard Haberman applied
> PDE's is nice, and if you want to see a nice new book with the cover
> showing solutions of PDE's plots which I am sure was made using
> Mathematica check Charles MacCluer's BVP and fourier expansions Dover
> book, but it does not contain any Mathematica code. Another book which
> uses a CAS system to solve PDE's is by David Betounes called PDE's for
> computational science, lots of examples and plots.
>
> Is there a short list somewhere which makes it clear what kind/class
> of PDE's Mathematica can solve and not solve directly using NDSolve or
> even DSolve? And why is it that NDSolve can solve an initial value PDE
> and not BVP? I wonder if NDSolve will be able to solve a BVP PDE in
> next version?
>
>
>>
Pratik Desai
>"James Gilmore" <james....@yale.edu> wrote in message
>news:dds9co$91c$1...@smc.vnet.net...
>
>
>
>>This is a classic mathematical physics BVP. You should approach this
>>problem
>>in Mathematica, as you would by hand: use separation of variables,
>>and then
>>a fourier expansion to satisfy the boundary conditions.
>>
>>
>
>Not if you want to use a numerical solvers such as NDSolve. That is
>the whole idea of using NDSolve.
>
>I know how to solve these by hand, and also by direct numerical
>approach, I've solved my of these before and more advanced ones when I
>took some courses at the Math dept at UC Berkeley one year ago, I was
>just playing around to see if NDSolve can solve BVP and get the same
>plots I got when I solved this problem by hand using sepration of
>variables.
>
>
>
>>There are many books
>>that explain how to do this.
>>
>>
>
>Yes, and my home library contains many fine books on PDE's. I like the
>Satnley Farlow book, and Mary Boas has excellent chapter on the
>subject, but a bit short on detailes. Also Richard Haberman applied
>PDE's is nice, and if you want to see a nice new book with the cover
>showing solutions of PDE's plots which I am sure was made using
>Mathematica check Charles MacCluer's BVP and fourier expansions Dover
>book, but it does not contain any Mathematica code. Another book which
>uses a CAS system to solve PDE's is by David Betounes called PDE's for
>computational science, lots of examples and plots.
>
>Is there a short list somewhere which makes it clear what kind/class
>of PDE's Mathematica can solve and not solve directly using NDSolve or
>even DSolve?
>
Would'nt that be awesome?
> And why is it that NDSolve can solve an initial value PDE
>and not BVP? I wonder if NDSolve will be able to solve a BVP PDE in
>next version?
>
>
Have you tried shooting methods, there must be a package doing just that
floating around
Best regards,
Pratik Desai
--
Pratik Desai
Graduate Student
UMBC
Department of Mechanical Engineering
Phone: 410 455 8134
James Gilmore
Thank you for your follow up, I will take a look at the books you suggested
when I have time. Apologies if you thought I suggested that you didn't know
how to solve Laplace like problems via hand. That was not intended.
There is a package on Math source that can be used to solve the Laplace
equation at least. I have it working on my machine (Mathematica 5.0.0.0 Win XP).
<http://library.wolfram.com/infocenter/MathSource/684>
The package needs development, but is a useful tool. I have solved the
Laplace equation with rectangular BC with this package. But beware, it has
trouble with corners (but I believe that can be fixed) and small numerical
errors can occur close to the boundaries.
--
James Gilmore
Graduate Student
Department of Physics
Yale University
New Haven, CT 06520 USA
"Nasser Abbasi" <n...@12000.org> wrote in message
news:dduroq$oih$1...@smc.vnet.net...
Mike Honeychurch
<n...@12000.org> wrote:
The NDSolve advanced documentation should be the first place to look for
this detailed information.
Cheers
Mike
Nasser Abbasi
>
> There is a package on Math source that can be used to solve the
> Laplace
> equation at least. I have it working on my machine (Mathematica
> 5.0.0.0 Win XP).
>
> <http://library.wolfram.com/infocenter/MathSource/684>
>
> The package needs development, but is a useful tool. I have solved
> the
> Laplace equation with rectangular BC with this package. But beware,
> it has
> trouble with corners (but I believe that can be fixed) and small
> numerical
> errors can occur close to the boundaries.
>
Thanks for the link to this package. It gives me errors just when
loading it. Here are the errors. This is on 5.2, windows XP.
http://12000.org/my_notes/LESolver/LESolver_errors.png
I also tried to execute the documentation notebook itself which
contains examples, parts of it worked, but many parts give many long
errors and it seems to have went in a loop, but may be I should have
waited a bit more.
Any way, since the package is dated 1996, I can understand. I am not
sure the original author is still at the email address shown there,
but I am sending them a link to the above now.
Nasser
Nasser Abbasi
news:ddscf6$ai1$1...@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
James Gilmore
Yes, I had similar bugs as well. These are just formatting errors in the
LESolver.m file. You need to go into the .m file and do
some manual editing. Most are spurious 'Returns' on the lines indicated.
I have been notified offline of another source:
<http://www.imtek.uni-freiburg.de/simulation/mathematica/IMSweb/>. The
server appears intermittent. I haven't tested the packages yet, but the
examples provided appear sound.
--
James Gilmore
Graduate Student
Department of Physics
Yale University
New Haven, CT 06520 USA
"Nasser Abbasi" <n...@12000.org> wrote in message
news:de45jl$qtq$1...@smc.vnet.net...
Peter Pein
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/
Nasser Abbasi
I've coded the solution to the 2-D steady state heat equation based on
a scheme called SOR (simultaneous over-relaxation), which is an
improved version of Gauss-Seidel scheme, which is in turn is an
improved version of the Jacobi method. These algorithms are described
in details in 'Numerical Methods for physics' by Garcia (nice book
btw, that was our text book for a course I took few years ago).
The Mathemetica implementation was straight forward and I have to say
that Mathematica 3D plots look really nice !
http://12000.org/my_notes/mma_matlab_control/e61/HTML/e61.htm
Nasser
Curtis Osterhoudt
on that file, and still haven't gotten it to work satifactorily. One of
the things that is NOT a formatting error is a function which is no
longer included past V. 3 of Mathematica: SingularValues. As far as I
can tell, this is replaced by the newer SingularValueList.
In any case, I'd love it if someone manages to get that package
working and posts it!
Cheers,
Curtis Osterhoudt
James Gilmore wrote:
>Hi Nasser,
>
>Yes, I had similar bugs as well. These are just formatting errors in the
>LESolver.m file. You need to go into the .m file and do
>some manual editing. Most are spurious 'Returns' on the lines indicated.
>
>I have been notified offline of another source:
><http://www.imtek.uni-freiburg.de/simulation/mathematica/IMSweb/>. The
>server appears intermittent. I haven't tested the packages yet, but the
>examples provided appear sound.
>
>
>
--
PGP Key ID: 0x235FDED1
Please avoid sending me Word or PowerPoint attachments.
http://www.gnu.org/philosophy/no-word-attachments.html
Curtis Osterhoudt
The problem of solving the heat equation with boundary conditions on a
plate has brought a lot of replies. "Mathematica Navigator", by Heikki
Ruskeepää, has a good section on such things in Chapter 23. I have the
first Edition, published in 1999, which is for version 3 of Mathematica.
I hear that there's been an updated version published since, and I'd
definitely get it if someone has had good experiences with it -- the 1st
Edition is an absolute gem, and is usually the book to which I turn
first if I don't know how to do something.
Section 23.1.1 "1D Parabolic and Hyperbolic Problems" starts with,
With NDSolve we can solve PDEs with one space and one time variable as the independent variables [...]
The command is not suitable for elliptic problems, or for problems with singularities. For elliptic problems
we present a finite difference method in Section 23.3.3.
In Section 23.3.3, Ruskeepää gives a quick function for solving 2D
elliptic problems, and provides a nice way of putting in the boundary
conditions, at least on a rectangular area, and the time-independent
source term, F(x,y). Here's (some of) the (transcribed) text for the
boundary conditions:
Consider the elliptic problem
u_xx + u_yy = F(x,y), x0 < x < x1, y0 < y < y1,
u(x, y1) = a2(x),
u(x0, y) = b1(y), u(x1, y) = b2(y),
u(x, y0) = a1(x).
Then, his function:
\!\(ellipticFD[{F_, \ a1_, \ a2_, \ b1_, \ b2_}, \ {x_, \ x0_, \ x1_, \
nx_}, \ {y_, \ y0_, \ y1_, \
ny_}] := \[IndentingNewLine]Module[{hx, \ hx2, \ xx, \ i, \ hy, \
hy2, \ yy, \ j, \ vars, \ u, \ eqs, \ c1, \ c2, \ c3, \ c4, \ sol1, \
sol2}, \ \[IndentingNewLine]\[IndentingNewLine]hx\ = \
N[\((x1\ - \ x0)\)/nx]; \ hx2\ = \ hx\^2; \
Do[xx[i]\ = \ x0\ + \ i\ hx, \ {i, \ 0, \
nx}]; \[IndentingNewLine]hy\ = \ N[\((y1\ - \ y0)\)/ny]; \
hy2\ = \ hy\^2; \
Do[yy[j]\ = \ y0\ + \ j\ hy, \ {j, \ 0, \
ny}]; \[IndentingNewLine]vars\ = \
Flatten[Array[u, \ {nx + 1, \ ny + 1}, \
0]]; \[IndentingNewLine]\[IndentingNewLine]eqs\ = \
Flatten[\[IndentingNewLine]\t\tTable[\[IndentingNewLine]\t\t\t\((u[
i\ + \ 1, \ j]\ - \ 2 u[i, j]\ + \
u[i - 1, \ j])\)/
hx2 + \[IndentingNewLine]\t\t\t\((u[i, \ j + 1]\ - \
2 u[i, \ j]\ + \ u[i, \ j\ - \ 1])\)/
hy2\ == \[IndentingNewLine]\t\t\tF /. {x -> xx[i], \
y -> yy[j]}, \ {i, \ nx\ - 1}, \ {j, \
ny\ - \ 1}]]; \[IndentingNewLine]\[IndentingNewLine]c1\ = \
Table[u[0\ \ , \ j\ ]\ == \ N[b1 /. y -> yy[j]], \ {j, \ 0, \
ny}]; \[IndentingNewLine]c2\ = \
Table[u[nx, j\ ]\ == \ N[b2 /. y -> yy[j]], \ {j, \ 0, \
ny}]; \[IndentingNewLine]c3\ = \
Table[u[i\ \ , 0\ ]\ == \ N[a1 /. x -> xx[i]], \ {i, \ 1,
nx\ - \ 1}]; \[IndentingNewLine]c4\ = \
Table[u[i\ \ , ny\ ]\ == \ N[a2 /. x -> xx[i]], \ {i, \ 1,
nx\ - \ 1}]; \[IndentingNewLine]\[IndentingNewLine]sol1\ = \
Solve[Join[c1, \ c2, \ c3, \ c4, \ eqs], \
vars]; \[IndentingNewLine]sol2\ = \
Table[u[i, j]\ /. \ sol1[\([1]\)], \ {i, \ 0, \ nx}, \ {j, \ 0, \
ny}]; \[IndentingNewLine]\[IndentingNewLine]ListInterpolation[
sol2, \ {{x0, \ x1}, \ {y0, \ y1}}]]\)
As his example, he says
We solve the elliptic problem
u_xx +u_yy = 0, 0<x<1, 0<y<1,
u(x,0) = 4x(1-x), u(0,y)=u(1,y)=u(x,1)=0.
Using 10 subintervals in each direction:
In[2]:= uappr = ellipticFD[{0, 4*x*(1 - x), 0, 0, 0}, {x, 0, 1, 10}, {y, 0, 1, 10}]
Out[2]=
InterpolatingFunction[]
In[3]:=
Plot3D[Evaluate[uappr[x, y]], {x, 0, 1}, {y, 0, 1}, ViewPoint -> {2., -2.4, 0.9}, AxesLabel -> {"x", "y", ""},
Ticks -> {{0, 1}, {0, 1}, {0, 1}}];
As is common in these problems, it's a bit tricky to match up boundary
conditions where corners don't have the same values for the conditions.
Still, I hope this works for someone!
Regards,
Curtis O.
Nasser Abbasi wrote:
--