(Uproszczone) zagadnienie Dirichleta dla prostokąta
 0⩽x≤ A, 0≤y≤B.
 Znaleźć funkcję u(x,y) spełniającą równanie Laplace’a 
 (((∂^2)u(x,y))/(∂(x^2))) + (((∂^2)u(x,y))/(∂(y^2))) = 0  (Δu(x,y)=0)
 oraz warunki brzegowe
 u(0,y) =ϕ(Subscript[, 1])(y),
 u(A,y) =ϕ(Subscript[, 2])(y),
 u(x,0) =ψ1(x),
 u(B,x) =ψ(Subscript[, 2])(x). 
 
Przyjmijmy najpierw, że ϕ(Subscript[, 1])(y)=ϕ(Subscript[, 2])(y)=0. Nazwijmy rozwiązanie dla tego przypadku vI(x,y). Potem można znaleźć rozwiązanie vII(x,y) dla ψ1(x)=ψ(Subscript[, 2])(x)=0. Pełne rozwiązanie będzie sumą rozwiązań vI i vII: v(x,y)=vI(x,y)+vII(x,y).
W tym notatniku skupimy się więc na poszukaniu vI(x,y).

A=2;

B=A/2;

(Subscript[ψ, 1])[x_]:=x*(A-x)

(Subscript[ψ, 2])[x_]:=-x^3*(A-x)

Plot[(Subscript[ψ, 1])[x],{x,0,A}]

Plot[(Subscript[ψ, 2])[x],{x,0,A}]

Próba użycia DSolve do rozwiązania równania
 (((∂^2)u(x,y))/(∂(x^2))) + (((∂^2)u(x,y))/(∂(y^2))) =0

rownanie=D[v[x,y],{x,2}]+D[v[x,y],{y,2}]==0;

warbrzeg={v[x,0]==(Subscript[ψ, 1])[x],v[x,B]==(Subscript[ψ, 2])[x],v[0,y]==0,v[A,y]==0};

sol11=DSolve[rownanie,v[x,y],{x,y}]

sol12=DSolve[Join[{rownanie},warbrzeg],v[x,y],{x,y}]

Od wersji 11 programu Mathematica dostajemy poprawne rozwiązanie w postaci nieskończonego szeregu !
Oczywiście w praktyce nie możemy sumować nieskończenie wielu składników i musimy uczynić sumę skończoną.

v0[x_,y_,kmax_]:=Sum[(-(1/((π^3) (k^3)))8 (-2+2 (((-1))^k)) Csch[(1/2) π k] Sin[(1/2) π x k] Sinh[(1/2) π (1-y) k]+
(1/((π^5) (k^5)))32 Csch[(1/2) π k] (24+6 (((-1))^k) (-4+(π^2) (k^2))) Sin[(1/2) π x k] Sinh[(1/2) π y k]),{k,1,kmax}]

N[v0[1/2,1/3,3]]

Plot3D[v0[x,y,1],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v0(x,y,1)",AxesLabel->{"x","y"}]

Plot3D[v0[x,y,20],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v0(x,y,20)",AxesLabel->{"x","y"}]

Próba użycia NDSolve do rozwiązania równania (((∂^2)u(x,y))/(∂(x^2))) + (((∂^2)u(x,y))/(∂(y^2))) =0

rownanie=D[v[x,y],{x,2}]+D[v[x,y],{y,2}]==0;

warbrzeg={v[x,0]==(Subscript[ψ, 1])[x],v[x,B]==(Subscript[ψ, 2])[x],v[0,y]==0,v[A,y]==0};

sol13=NDSolve[Join[{rownanie},warbrzeg],v,{x,0,A},{y,0,B}]

v1[x_,y_]=v[x,y]/.sol13[[1]]

Plot3D[v1[x,y],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v1(x,y)",AxesLabel->{"x","y"}]

Od wersji 11 programu Mathematica można bezpośrednio rozwiązać równanie Laplace’a przy pomocy NDSolve!

W wersjach poprzednich było to niemożliwe, więc należało zbudować  rozwiązanie z wykorzystaniem metody relaksacji (poprzez wprowadzenie pomocniczej zależności od czasu), gdy równanie przyjmuje postać:
((∂v(t,x,y))/(∂t)) - κ( (((∂^2)v(t,x,y))/(∂(x^2))) + (((∂^2)v(t,x,y))/(∂(y^2))) ) = 0.
Dla wystarczająco dużych czasów
funkcja v(t,x,y) przestaje się zmieniać z czasem i zachodzi u(x,y)=v((Subscript[t, max]),x,y).

rownanie=D[v[t,x,y],t]-κ*(D[v[t,x,y],{x,2}]+D[v[t,x,y],{y,2}])==0;

warbrzeg={v[t,x,0]==(Subscript[ψ, 1])[x],v[t,x,B]==(Subscript[ψ, 2])[x],v[t,0,y]==0,v[t,A,y]==0};

warpocz=v[0,x,y]==(Subscript[ψ, 1])[x]+((Subscript[ψ, 2])[x]-(Subscript[ψ, 1])[x])/B*y;

tmax=2;κ=1;

sol=NDSolve[Join[{rownanie},warbrzeg,{warpocz}],v,{t,0,tmax},{x,0,A},{y,0,B}]

v[t_,x_,y_]=v[t,x,y]/.sol[[1]]

Plot3D[v[0,x,y],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v(0,x,y)",AxesLabel->{"x","y"}]

Plot[v[0,A/2,y],{y,0,B},PlotRange->All,AxesLabel->{"y","v(0,A/2,y)"}]

Plot3D[v[tmax/2,x,y],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v(tmax/2,x,y)",AxesLabel->{"x","y"}]

Plot[v[tmax/2,A/2,y],{y,0,B},PlotRange->All,AxesLabel->{"y","v(tmax/2,A/2,y)"}]

Plot3D[v[tmax,x,y],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v(tmax,x,y)",AxesLabel->{"x","y"}]

Plot[v[tmax,A/2,y],{y,0,B},PlotRange->All,AxesLabel->{"y","v(tmax,A/2,y)"}]

Plot3D[v[tmax,x,y]-v[tmax-1/10,x,y],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"v(tmax,x,y)-v(tmax-1/10,x,y)",AxesLabel->{"x","y"}]

Rozwiązanie analityczne z wykorzystaniem szeregów Fouriera jest równoważne temu, co dostaliśmy wcześniej przy użyciu funkcji DSolve

u[x_,y_,n_]:=(a[n]*Sinh[n*Pi/A*(B-y)]+b[n]*Sinh[n*Pi/A*y])*Sin[n*Pi/A*x];

Simplify[D[u[x,y,n],{x,2}]+D[u[x,y,n],{y,2}]]

u[0,y,n]

Simplify[u[A,y,n],Element[n,Integers]]

a[n_]=2/(A*Sinh[n*Pi*B/A])*Integrate[(Subscript[ψ, 1])[x]*Sin[n*Pi/A*x],{x,0,A}]

b[n_]=2/(A*Sinh[n*Pi*B/A])*Integrate[(Subscript[ψ, 2])[x]*Sin[n*Pi/A*x],{x,0,A}]

U[x_,y_,nmax_]:=Sum[u[x,y,n],{n,1,nmax}]

U[x,y,3]

Plot3D[Evaluate[U[x,y,3]],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"nmax=3",AxesLabel->{"x","y"}]

Plot3D[Evaluate[U[x,y,5]],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"nmax=5",AxesLabel->{"x","y"}]

Plot3D[Evaluate[U[x,y,7]],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"nmax=7",AxesLabel->{"x","y"}]

Plot3D[Evaluate[U[x,y,30]],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"nmax=30",AxesLabel->{"x","y"},PlotLabel->"V(x,y)"]

Plot3D[Evaluate[U[x,y,51]-U[x,y,50]],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"U(nmax=51)-U(nmax=50)",AxesLabel->{"x","y"}]

Rozwiązanie numeryczne w oparciu o schemat Jacobiego

(* Siatka prostokątna *)

Nx=40;Δx=A/Nx;

Ny=Nx/2;Δy=B/Ny;

Do[x[i]=i*Δx,{i,0,Nx}]

Do[y[j]=j*Δy,{j,0,Ny}]

(* Warunki brzegowe *)

Do[US[0,j]=0,{j,0,Ny}]

Do[US[Nx,j]=0,{j,0,Ny}]

Do[US[i,0]=(Subscript[ψ, 1])[i*Δx],{i,0,Nx}]

Do[US[i,Ny]=(Subscript[ψ, 2])[i*Δx],{i,0,Nx}]

(* Sytuacja początkowa wewnątrz prostokąta *)

Do[US[i,j]=1,{i,1,Nx-1},{j,1,Ny-1}]

dane0=Table[{x[i],y[j],US[i,j]},{i,0,Nx},{j,0,Ny}];

(* Iteracje *)

Niter=100;

Do[USSTARY[i,j]=US[i,j],{i,0,Nx},{j,0,Ny}]

Do[{Do[USNOWY[i,j]=(USSTARY[i,j-1]+USSTARY[i,j+1]+USSTARY[i-1,j]+USSTARY[i+1,j])/4,{i,1,Nx-1},{j,1,Ny-1}];Do[USSTARY[i,j]=USNOWY[i,j],{i,1,Nx-1},{j,1,Ny-1}]},{iter,1,Niter}]

dane100=Table[{x[i],y[j],USNOWY[i,j]},{i,0,Nx},{j,0,Ny}];

ListPlot3D[{dane100},PlotLabel->"US100",AxesLabel->{"x","y"}]

(* Dodatkowe iteracje *)

Niter=100;

Do[{Do[USNOWY[i,j]=(USSTARY[i,j-1]+USSTARY[i,j+1]+USSTARY[i-1,j]+USSTARY[i+1,j])/4,{i,1,Nx-1},{j,1,Ny-1}];Do[USSTARY[i,j]=USNOWY[i,j],{i,1,Nx-1},{j,1,Ny-1}]},{iter,1,Niter}]

dane200=Table[{x[i],y[j],USNOWY[i,j]},{i,0,Nx},{j,0,Ny}];

ListPlot3D[{dane200},PlotLabel->"US200",AxesLabel->{"x","y"}]

(* I jeszcze 100 nowych iteracji *)

Niter=100;

Do[{Do[USNOWY[i,j]=(USSTARY[i,j-1]+USSTARY[i,j+1]+USSTARY[i-1,j]+USSTARY[i+1,j])/4,{i,1,Nx-1},{j,1,Ny-1}];Do[USSTARY[i,j]=USNOWY[i,j],{i,1,Nx-1},{j,1,Ny-1}]},{iter,1,Niter}]

dane300=Table[{x[i],y[j],USNOWY[i,j]},{i,0,Nx},{j,0,Ny}];

ListPlot3D[{dane300},PlotLabel->"US300",AxesLabel->{"x","y"}]

(* I jeszcze 200 iteracji *)

Niter=200;

Do[{Do[USNOWY[i,j]=(USSTARY[i,j-1]+USSTARY[i,j+1]+USSTARY[i-1,j]+USSTARY[i+1,j])/4,{i,1,Nx-1},{j,1,Ny-1}];Do[USSTARY[i,j]=USNOWY[i,j],{i,1,Nx-1},{j,1,Ny-1}]},{iter,1,Niter}]

dane500=Table[{x[i],y[j],USNOWY[i,j]},{i,0,Nx},{j,0,Ny}];

ListPlot3D[{dane500},PlotLabel->"US500",AxesLabel->{"x","y"}]

danecomp=Table[{x[i],y[j],-USNOWY[i,j]+U[i*Δx,j*Δy,50]},{i,0,Nx},{j,0,Ny}];

ListPlot3D[{danecomp},PlotLabel->"Uana(50)-Unum(500)",AxesLabel->{"x","y"}]

Plot3D[Evaluate[U[x,y,50]-v[tmax,x,y]],{x,0,A},{y,0,B},PlotRange->All,PlotLabel->"U(nmax=50)-v(tmax,x,y)",AxesLabel->{"x","y"}]