Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Could some one please help

• To: mathgroup at smc.vnet.net
• Subject: [mg99945] Re: [mg99924] Could some one please help
• 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[
    Flatten[MapThread[
      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

• References:
  • Could some one please help
    • From: gtw <graceneetw@google.com>