       • To: mathgroup at smc.vnet.net
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 20 May 2009 04:56:40 -0400 (EDT)
• References: <200905191046.GAA06635@smc.vnet.net>

```gtw wrote:
> I am new in numerical.  The following is
> to try to solve D[u[x,y],x,x]+D[u[x,y],y,y]=1-u[x,y]
> The line eqns=... is wrong (I think)
> Could someone tell me how to fix?
> Besides, how to put Neumann boundary condition
> Thanks
>
> num = 80; hStep = N[1/(num + 1)];
> vars = Table[u[i, j], {i, num}, {j, num}];
> Short[vars, 3]
> kern = {{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}/hStep^2;
> lap = ListCorrelate[kern, vars, {2, 2}, 0];
> (* equation for solution on grid*)Short[
> eqns = Thread[Map[Flatten, lap == 1 - u[i, j]]], 20]
> {vec, mat} = CoefficientArrays[eqns, Flatten[vars]]
> sol = LinearSolve[mat, -vec];
> sol = Partition[sol, num];
> ListContourPlot[sol]
> ListPlot3D[sol]

The variant below should set up the equations correctly.

(*equation for solution on grid*)
eqns = Flatten[
MapThread[Equal, {lap, Table[1 - u[i, j], {i, num}, {j, num}]},
2]];

Also, I suspect you want:

kern = {{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}*hStep^2;

(Note use of multiplication instead of division.)

To attain Neumann or other boundary conditions you would need to remove
the boundary equations from the set above. Then replace them with
suitable approximations to the boundary conditions you wish to impose.
But you need not actually remove anything because you can simply set up
so that your lap only considers the inside rectangle.

Example:

lap = ListCorrelate[kern, vars];

bdy = Join[Table[u[1, j] == j, {j, num}],
Table[u[num, j] == j^2, {j, num}],
Table[u[j, 1] == Sin[j], {j, 2, num - 1}],
Table[u[j, num] == j^3, {j, 2, num - 1}]];

eqns = Join[
Equal, {lap,
Table[1 - u[i, j], {i, 2, num - 1}, {j, 2, num - 1}]}, 2]],
bdy];

Of course one could instead use derivative approximations via finite
differences for the boundary, or mix linear combinations of boundary
values plus boundary derivative values.

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: Re: Solving for Stiffness using a plot
• Next by Date: Re: Customise a Mathematica 7.0 Palette for braces and square brackets