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: [mg59812] Re: [mg59645] How to specify boundary conditions on all 4 sides of a plate for a steady state heat equation (PDE) using NDSolve? (Laplace equation)
- From: Curtis Osterhoudt <gardyloo at mail.wsu.edu>
- Date: Mon, 22 Aug 2005 02:48:39 -0400 (EDT)
- References: <200508151050.GAA25179@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi, all,
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:
>
>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
>
>
>
>
>
>
>
>
>
>
--
PGP Key ID: 0x235FDED1
Please avoid sending me Word or PowerPoint attachments.
http://www.gnu.org/philosophy/no-word-attachments.html
- References:
- How to specify boundary conditions on all 4 sides of a plate for a steady state heat equation (PDE) using NDSolve? (Laplace equation)
- From: "Nasser Abbasi" <nma@12000.org>
- How to specify boundary conditions on all 4 sides of a plate for a steady state heat equation (PDE) using NDSolve? (Laplace equation)