# Solving Laplace equation

I tried to solve the Laplace equation by two different methods, I made two 3D plots, but they didn't look the same. Is anything wrong? what can I do? I tried to use Fourier series at first and some Mathematica code

Clear["Global*"]
(*rectangle 0<x<L,0<y<H*)
L = 4;
H = 3;

(*boundary conditions*)
f1[x_] := If[L > x > 0, 23]
f2[x_] := If[L > x > 0, 1]
g1[y_] := If[H > y > 0, -2]
g2[y_] := If[H > y > 0, 31]

(*coefficient A(n)*)
A[n_] := (2/(L*Sinh[n*Pi*H/L]))*Integrate[f2[x]*Sin[(n*Pi*x)/L], {x, 0, L}]
A[n]

(*coefficient B(n)*)
B[n_] := (2/(L*Sinh[n*Pi*H/L]))*Integrate[f1[x]*Sin[(n*Pi*x)/L], {x, 0, L}]
B[n]

(*coefficient C(n)*)
C1[n_] := (2/(H*Sinh[n*Pi*L/H]))*Integrate[g2[y]*Sin[(n*Pi*y)/L], {y, 0, H}]
C1[n]

(*coefficient D(n)*)
D1[n_] := (2/(H*Sinh[n*Pi*L/H]))*Integrate[g1[y]*Sin[(n*Pi*y)/L], {y, 0, H}]
D1[n]

(*Partial solutions*)
w1[x_, y_, N_] := Sum[A[n]*Sinh[n*Pi*y/L]*Sin[n*Pi*x/L], {n, 1, N}]
w2[x_, y_, N_] := Sum[B[n]*Sinh[n*Pi*y/L]*Sin[n*Pi*x/L], {n, 1, N}]
w3[x_, y_, N_] := Sum[C1[n]*Sinh[n*Pi*x/H]*Sin[n*Pi*y/H], {n, 1, N}]
w4[x_, y_, N_] := Sum[D1[n]*Sinh[n*Pi*x/H]*Sin[n*Pi*y/H], {n, 1, N}]
w[x_, y_, N_] := w1[x, y, N] + w2[x, y, N] + w3[x, y, N] + w4[x, y, N]

(*General solution*)
sol2 = w[x, y, 30]

Plot3D[sol2, {x, 0, L}, {y, 0, H}, Mesh -> Automatic,
MeshFunctions -> {#3 &}, PlotTheme -> Automatic,
AxesLabel -> {"x", "y", "w"}, PlotLabel -> "Laplace equation"]

leqn = Laplacian[u[x, y], {x, y}] == 0;
bc = {u[x, 0] == 23, u[x, H] == 1, u[0, y] == -2, u[L, y] == 31};
sol = FullSimplify[u[x, y] /. DSolve[{leqn, bc}, u[x, y], {x, y}][[1]]]
asol = sol /. \[Infinity] -> 30 // Activate
Plot3D[asol, {x, 0, L}, {y, 0, H}, PlotRange -> All]

• A question after a quick look - are your 'boundary conditions' (f1, f2, g1, g2) in the first section correct? They seem to be set to constants across the whole rectangle and are independent of the alternate variable, I.e. f1 and f2 do not depend on y and do not appear to be the boundary. As I said this is just a quick look while waiting for a plane.
– jmm
Commented Apr 6, 2023 at 19:07
• I tried some constants for simplicity, I tried and some functions too Commented Apr 6, 2023 at 19:08
• If I understand this correctly, then the boundary conditions at the vertices are incompatible. Commented Apr 6, 2023 at 19:13
• Ι have wrong somewhere, I couldn't find it Commented Apr 6, 2023 at 19:16

There is an error or typo in the scaling of the arguement in y-direction with H instead of L. Your If-Functions don't produce a value outside the True-Interval

If[1 > 2, 0]
??


Both coefficient definitions don't match boundary values in the second method with NDSolve

C1[n_] := (2/(H*Sinh[n*Pi*L/H]))*Integrate[g2[y]*Sin[(n*Pi*y)/L], {y, 0, H}]

(*coefficient D(n)*)

D1[n_] := (2/(H*Sinh[n*Pi*L/H]))*Integrate[g1[y]*Sin[(n*Pi*y)/L], {y, 0, H}]
`

Probably it would be worth to try a { L -> H } because the boundary functions in the second are constant. Or the other way around.

Another bad idea is to use single captital letters as Symbols, even in Patterns, as you are aware of. You didnt make use of C,D but there is a N.

The definition with Pattern N yields a silly error message if you try to Simplify the expression.