NotebookDirectory[]

Equilibrium temperature distribution - 2D problem

The problem

The objective is to solve the PDE 

(((∂^2)u)/(∂(x^2)))(x,y)+(((∂^2)u)/(∂(y^2)))(x,y)=0  on  {(x,y) ∈  | 0≤x ≤K, 0 ≤y≤L},

subject to the conditions

u(x,0) = (Subscript[f, 1])(x), u(x,L)=(Subscript[f, 2])(x)   (call these BCx) 

u(0,y) = (Subscript[g, 1])(y), u(K,y)=(Subscript[g, 2])(y)  (call these BCy)

The trick is to split this problem into two problems

Problem 1

The objective is to solve the PDE 

(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))=0  on  0≤x ≤K, 0 ≤y≤L,

subject to the conditions

u(x,0) = (Subscript[f, 1])(x), u(x,L)=(Subscript[f, 2])(x)  (call these BCx) 

u(0,y) = 0, u(K,y)=0  (call these DBCy)

Step 1. First ignore BCx  and solve

(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))=0  on  0≤x ≤K, 0 ≤y≤L,

subject to the conditions

u(0,y) = 0, u(K,y)=0  (call these DBCy)

Using SofV method we find few solutions of this problem:

Sin[((n Pi)/K)x] ((Sinh[((n Pi)/K)y])/(Sinh[((n Pi)/K)L]))

Sin[((n Pi)/K)x] ((Sinh[((n Pi)/K)(L-y)])/(Sinh[((n Pi)/K)L]))

Test these solutions:

(D[#,{x,2}]+D[#,{y,2}])&[Sin[((n Pi)/K)x](( Sinh[((n Pi)/K)y])/( Sinh[((n Pi)/K)L]))]

FullSimplify[(Sin[((n Pi)/K)x] (( Sinh[((n Pi)/K)y])/( Sinh[((n Pi)/K)L])))/.{x->{0,K}}]

We need to tell Mathematica that n is an integer.

FullSimplify[(Sin[((n Pi)/K)x] (( Sinh[((n Pi)/K)y])/( Sinh[((n Pi)/K)L])))/.{x->{0,K}},n∈Integers]

(D[#,{x,2}]+D[#,{y,2}])&[Sin[((n Pi)/K)x] ((Sinh[((n Pi)/K)(L-y)])/( Sinh[((n Pi)/K)L]))]

FullSimplify[(Sin[((n Pi)/K)x](( Sinh[((n Pi)/K)(L-y)])/(Sinh[((n Pi)/K)L])))/.{x->{0,K}}]

FullSimplify[(Sin[((n Pi)/K)x](( Sinh[((n Pi)/K)(L-y)])/(Sinh[((n Pi)/K)L])))/.{x->{0,K}},n∈Integers]

Step 2. Now that we have few solutions we form many solutions.  This is the FFM idea, which is commonly known as the superposition principle:

(∑_(n=1)^nn) (Subscript[a, n ])Sin[((n Pi)/K)x] ((Sinh[((n Pi)/K)y])/(Sinh[((n Pi)/K)L]))+(∑_(n=1)^nn) (Subscript[b, n ])Sin[((n Pi)/K)x](( Sinh[((n Pi)/K)(L-y)])/(Sinh[((n Pi)/K)L]))

Next we choose (Subscript[a, n]) and (Subscript[b, n]) such that the above function satisfies BCx conditions. First substitute y=0. This leads to the formula for (Subscript[b, n]).

(Subscript[f, 1])(x) =(∑_(n=1)^nn) (Subscript[b, n ])Sin[((n Pi)/K)x]

(Subscript[f, 1])(x)Sin[((j Pi)/K)x] =(∑_(n=1)^nn) (Subscript[b, n ])Sin[((n Pi)/K)x]Sin[((j Pi)/K)x]

Clear[K]; FullSimplify[Integrate[Sin[((n Pi)/K)x]Sin[((j Pi)/K)x],{x,0,K}],And[n∈Integers,j∈Integers,Or[j>n,j<n]]]

Clear[K]; FullSimplify[Integrate[Sin[((n Pi)/K)x]Sin[((n Pi)/K)x],{x,0,K}],And[n∈Integers]]

(∫_0^K)(Subscript[f, 1])(x)Sin[((j Pi)/K)x]dx =(∑_(n=1)^nn) (Subscript[b, n ])(∫_0^K)Sin[((n Pi)/K)x]Sin[((j Pi)/K)x]dx

(∫_0^K)(Subscript[f, 1])(x)Sin[((j Pi)/K)x]dx = (Subscript[b, n ])(∫_0^K)Sin[((j Pi)/K)x]Sin[((j Pi)/K)x]dx

(∫_0^K)(Subscript[f, 1])(x)Sin[((j Pi)/K)x]dx = (Subscript[b, j ])*(K/2)

Then substitute y=L. This leads to the formula for (Subscript[a, n]).

Problem 2

The objective is to solve the PDE 

(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))=0  on  0≤x ≤K, 0 ≤y≤L,

subject to the conditions

u(x,0) = 0, u(x,L)=0  (call these DBCx) 

u(0,y) = (Subscript[g, 1])(y), u(K,y)=(Subscript[g, 2])(y)  (call these BCy)

Step 1. First ignore BCy  and solve

(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))=0  on  0≤x ≤K, 0 ≤y≤L,

subject to the conditions

u(x,0) = 0, u(x,L)=0  (call these DBCx)

Using SofV method we find few solutions of this problem:

Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)x])/(Sinh[((n Pi)/L)K]))

Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)(K-x)])/(Sinh[((n Pi)/L)K]))

Test these solutions:

(D[#,{x,2}]+D[#,{y,2}])&[Sin[((n Pi)/L)y](( Sinh[((n Pi)/L)x])/(Sinh[((n Pi)/L)K]))]

FullSimplify[(Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)x])/(Sinh[((n Pi)/L)K])))/.{y->{0,L}},n∈Integers]

(D[#,{x,2}]+D[#,{y,2}])&[Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)(K-x)])/(Sinh[((n Pi)/L)K]))]

FullSimplify[(Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)(K-x)])/(Sinh[((n Pi)/L)K])))/.{y->{0,L}},n∈Integers]

Step 2. Now that we have few solutions we form many solutions.  This is the FFM idea, which is commonly known as the superposition principle:

(∑_(n=1)^nn) (Subscript[c, n ])Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)x])/(Sinh[((n Pi)/L)K]))+(∑_(n=1)^nn) (Subscript[d, n ])Sin[((n Pi)/L)y] ((Sinh[((n Pi)/L)(K-x)])/(Sinh[((n Pi)/L)K]))

Now we choose (Subscript[c, n]) and (Subscript[d, n]) such that the above function satisfies BCy conditions. First substitute x=0. This leads to the formula for (Subscript[d, n]). Then substitute x=K. This leads to the formula for (Subscript[c, n]).

A symbolic implementation

Here are the given quantities

Clear[lK,lL,f1,f2,g1,g2,nn];

nn=15; 

lK = 1; lL = 1;

 f1[x_]=4(x^2)(1-x); 

g2[y_]=4y(((1-y))^2); 

f2[x_]=0; 

g1[y_]=0;

Clear[aa];

aa[n_]=(2/(lK ))Integrate[f2[x]Sin[((n Pi)/lK)x],{x,0,lK}]

Clear[bb];

bb[n_]=(2/(lK ))Integrate[f1[x]Sin[((n Pi)/lK)x],{x,0,lK}]

Clear[cc];

cc[n_]=(2/(lL ))Integrate[g2[y]Sin[((n Pi)/lL)y],{y,0,lL}]

Clear[dd];

dd[n_]=(2/(lL ))Integrate[g1[y]Sin[((n Pi)/lL)y],{y,0,lL}]

The solution is

Clear[uu]; 

uu[x_,y_]=(∑_(n=1)^nn) aa[n]Sin[((n Pi)/lK)x] ((Sinh[((n Pi)/lK)y])/(Sinh[((n Pi)/lK)lL]))+(∑_(n=1)^nn) bb[n]Sin[((n Pi)/lK)x] ((Sinh[((n Pi)/lK)(lL-y)])/(Sinh[((n Pi)/lK)lL]))+ (∑_(n=1)^nn) cc[n]Sin[((n Pi)/lL)y] ((Sinh[((n Pi)/lL)x])/(Sinh[((n Pi)/lL)lK]))+(∑_(n=1)^nn) dd[n] Sin[((n Pi)/lL)y] ((Sinh[((n Pi)/lL)(lK-x)])/(Sinh[((n Pi)/lL)lK]));

How good is our approximation for the function f1[x] in the boundary conditions? Here is a visual answer.

Plot[{f1[x],uu[x,0]},{x,0,1},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.005]}}]

Plot[uu[x,lL],{x,0,1}]

Plot[uu[0,x],{x,0,1}]

Plot[{g2[y],uu[lK,y]},{y,0,1},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.005]}}]

Plot3D[N[uu[x,y]],{x,0,1},{y,0,1},Mesh->Automatic]

DensityPlot[N[uu[x,y]],{x,0,1},{y,0,1},Frame->True,PlotRange->{{0,1},{0,1}},ColorFunction->(RGBColor[1,1-#,1-#]&)]

A numerical implementation

Here are the given quantities

Clear[lK,lL,f1,f2,g1,g2,nn];

nn=45; 

lK = Pi; lL = Pi;

 f1[x_]=(((Sin[x]))^(1/10));g2[y_]=0; f2[x_]=0; g1[y_]=0;

Plot[f1[x],{x,0,lK},PlotRange->{0,1.2}]

Clear[aa];

aa[n_]=(2/lK)Integrate[f2[x]Sin[((n Pi)/lK)x],{x,0,lK}]

Clear[bbl];

bbl=Chop[Table[(2/lK)NIntegrate[f1[x]Sin[((n Pi)/lK)x],{x,0,lK},
Method->{Automatic},MaxRecursion->200,AccuracyGoal->12,PrecisionGoal->16],{n,1,nn}]]

Clear[cc];

cc[n_]=(2/(lL Sinh[((n Pi)/lL)lK]))Integrate[g2[y]Sin[((n Pi)/lL)y],{y,0,lL}]

Clear[dd];

dd[n_]=(2/(lL ))Integrate[g1[y]Sin[((n Pi)/lL)y],{y,0,lL}]

The solution is

Clear[uu];

uu[x_,y_]=(∑_(n=1)^nn) aa[n]Sin[((n Pi)/lK)x] ((Sinh[((n Pi)/lK)y])/(Sinh[((n Pi)/lK)lL]))+(∑_(n=1)^nn) bbl[[n]]Sin[((n Pi)/lK)x] ((Sinh[((n Pi)/lK)(lL-y)])/(Sinh[((n Pi)/lK)lL]))+ (∑_(n=1)^nn) cc[n]Sin[((n Pi)/lL)y] ((Sinh[((n Pi)/lL)x])/(Sinh[((n Pi)/lL)lK]))+(∑_(n=1)^nn) dd[n] Sin[((n Pi)/lL)y] ((Sinh[((n Pi)/lL)(lK-x)])/(Sinh[((n Pi)/lL)lK]));

Plot[{f1[x],uu[x,0]},{x,0,lK},PlotRange->{0,1.2}]

Plot[uu[x,lL],{x,0,1}]

Plot[uu[0,x],{x,0,1}]

Plot[uu[lK,x],{x,0,1}]

Plot3D[N[uu[x,y]],{x,0,lK},{y,0,lL}]

DensityPlot[N[uu[x,y]],{x,0,lK},{y,0,lL},Frame->True,PlotRange->{{0,lK},{0,lL}},ColorFunction->(RGBColor[1,1-#,1-#]&)]

A symbolic implementation with a problem

Here are the given quantities

Clear[lK,lL,f1,f2,g1,g2,nn];

nn=20; 

lK = 1; lL = 1;

 f1[x_]=3(x^2)+1; 

g2[y_]=4-8y (y-1); 

 f2[x_]=4+8 x (x-1); 

g1[y_]=1+3(y^2);

Clear[aa];

aa[n_]=(2/(lK ))Integrate[f2[x]Sin[((n Pi)/lK)x],{x,0,lK}]

Clear[bb];

bb[n_]=(2/(lK ))Integrate[f1[x]Sin[((n Pi)/lK)x],{x,0,lK}]

Clear[cc];

cc[n_]=(2/(lL ))Integrate[g2[y]Sin[((n Pi)/lL)y],{y,0,lL}]

Clear[dd];

dd[n_]=(2/(lL ))Integrate[g1[y]Sin[((n Pi)/lL)y],{y,0,lL}]

The solution is

Clear[uu]; 

uu[x_,y_]=(∑_(n=1)^nn) aa[n]Sin[((n Pi)/lK)x] ((Sinh[((n Pi)/lK)y])/(Sinh[((n Pi)/lK)lL]))+(∑_(n=1)^nn) bb[n]Sin[((n Pi)/lK)x] ((Sinh[((n Pi)/lK)(lL-y)])/(Sinh[((n Pi)/lK)lL]))+ (∑_(n=1)^nn) cc[n]Sin[((n Pi)/lL)y] ((Sinh[((n Pi)/lL)x])/(Sinh[((n Pi)/lL)lK]))+(∑_(n=1)^nn) dd[n] Sin[((n Pi)/lL)y] ((Sinh[((n Pi)/lL)(lK-x)])/(Sinh[((n Pi)/lL)lK]));

How good are the approximations? Here are visual answers:

Plot[{f1[x],uu[x,0]},{x,0,1},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.005]}},PlotRange->All]

Plot[{g2[y],uu[lK,y]},{y,0,1},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.005]}},PlotRange->All]

Plot[{f2[x],uu[x,lL]},{x,0,1},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.005]}},PlotRange->All]

Plot[{g1[y],uu[0,y]},{y,0,1},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.005]}},PlotRange->All]

Plot3D[N[uu[x,y]],{x,0,1},{y,0,1},Mesh->Automatic,PlotRange->{0,6.5}]

DensityPlot[N[uu[x,y]],{x,0,1},{y,0,1},Frame->True,PlotRange->{{0,1},{0,1}},ColorFunction->(RGBColor[1,1-#,1-#]&)]