Hi,

I solve laplace equation in a square. All the lines of my code is okay.

Please just look to the last part of my code titled procedure:

When I run my code without (last funciton f #f := (x,y) -> 0;) see last lines to find "f". It's runing, there is no problem. But when I put add f, there is an error. Many think for any help.

Procedure

Using the previous suty in section stencil we can write the procedure to solve the Laplace equation in [0,1]*[0,1] with the boundary condition Neumann conditions on the vertical boundary and Dirichlet boundary condition on the horizontal baoundary. In our study we will use the same stepsize h in x and y direction.

PoissonSolve:=proc(N,_f)

local Z,i,h,y,x,sys,w,f,sol,j,u,Data;

# define basic grid parameters

Z := i -> (1/(N+1))*i;

x[0] = Z(0),x[N+1] = Z(N+1),y[0] = Z(0),y[N+1] = Z(N+1);

h := evalf(Z(1)-Z(0));

# Fix the boundary data and the source matrix

for i from 0 to N+1 do:

# Neumann boundary condition

u[N+1,i] := u[N,i] ;

u[0,i] := u[1,i];

# Dirichlet boundary condition

u[i,0] := 0;

u[i,N+1] := 0;

od:

f := Array(0..N+1,0..N+1,[seq([seq(evalf(_f(Z(i),Z(j))),i=0..N+1)],j=0..N+1)],datatype=float);

# Write down the system of equations to solve and solve them

sys := [seq(seq(Stencil(h,i,j,u,f),i=1..N),j=1..N)];

w := [seq(seq(u[i,j],i=1..N),j=1..N)];

sol := LinearSolve(GenerateMatrix(sys,w));

# parse the solution vector sol back into "matrix" form

for i from 1 to N do:

for j from 1 to N do:

u[i,j] := sol[(j-1)*N+i]:

od:

od:

# generate a 3D plot of the solution using the surfdata command

Data := [seq([seq([Z(i),Z(j),u[i,j]],i=0..N+1)],j=0..N+1)]:

surfdata(Data,axes=boxed,labels=[`x`,`y`,`u(x,y)`],shading=zhue,style=patchcontour);

end proc:

Here is an example of the output when the source function is set to zero

f(x, y) = 0

; i.e., when reduces down to Laplace's equation:

#f := (x,y) -> 0;

#PoissonSolve(10,f);

Question8.mw