# Simplifying the final coefficient expression of the solution to a Laplacian

I have been trying to solve a boundary value problem analytically which involves a three-dimensional temperature Laplacian over a parallelepiped. In the final step of my solution, using the two non-homogeneous $$z$$-boundary conditions, I calculate the two unknown Fourier coefficients $$C_1,C_2$$. The mathematica code is as follows:

T[x_, y_, z_] = (C1*E^(γ z) + C2*E^(-γ z))*Sin[(α x/L) + β]*Sin[(δ y/l) + θ] + Ta;

tc[x_, y_] = E^(-bc*y/l)*{tci + (bc/l)*Integrate[E^(bc*s/l)*T[x, s, 0], {s, 0, y}]};

bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);

ortheq1 = Integrate[(bc1[] - bc1[])*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}, Assumptions -> {C1 > 0, C2 > 0, L > 0, l > 0, α > 0, β > 0, γ > 0, δ > 0, θ > 0, bc > 0, pc > 0, Ta > 0, tci > 0}] == 0;

th[x_, y_] = E^(-bh*x/L)*{thi + (bh/L)*Integrate[E^(bh*s/L)*T[s, y, w], {s, 0, x}]};

bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);

ortheq2 = Integrate[(bc2[] - bc2[])*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}, Assumptions -> {C1 > 0, C2 > 0, L > 0, l > 0, α > 0, β > 0, γ > 0, δ > 0, θ > 0, bh > 0, ph > 0, Ta > 0, thi > 0}] == 0;

soln = Solve[{ortheq1, ortheq2}, {C1, C2}];

CC1 = C1 /. soln[[1, 1]];
CC2 = C2 /. soln[[1, 2]];



The coefficients CC1, CC2 I get from this code are extremely complex and lengthy (I must mention here that they are correct as I have verified my series solution against an FEM approach), which makes reporting them in a thesis or a scientific communication troublesome.

I will appreciate if someone can help me simplify the resulting C1,C2 expressions. I have already tried the inbuilt Simplify command to not so favorable results.

Some context and possibly helpful information

The function I am trying to determine is of the form:

$$T(x,y,z)=\sum_{n,m=0}^{\infty}(C_1 e^{\gamma z}+C_2 e^{-\gamma z})\sin\bigg(\frac{\alpha_n x}{L}+\beta_n\bigg)\sin\bigg(\frac{\delta_m y}{l}+\theta_m\bigg)+T_a \tag 1$$

The two $$z$$ bc(s) are:

$$\frac{\partial T(x,y,0)}{\partial z}=p_c (T(x,y,0)-t_c) \tag 2$$ $$\frac{\partial T(x,y,w)}{\partial z}=p_h (t_h-T(x,y,w))\tag 3$$

I have defined $$t_c, t_h$$ in my code and am not repeating here. If someone would have solved this problem using a pen and paper approach, he/she would have substituted $$(1)$$ in $$(2), (3)$$ and multiplied throughout by $$\int_0^L \sin\bigg(\frac{\alpha_k x}{L}+\beta_k\bigg)\mathrm{d}x$$ and $$\int_0^l \sin\bigg(\frac{\delta_j y}{l}+\theta_j\bigg)\mathrm{d}y$$ and used their orthogonality to remove the summations. In this procedure, he/she might have used the following relations:

$$u=\int_0^L \sin\bigg(\frac{\alpha_n x}{L}+\beta_n\bigg)\sin\bigg(\frac{\alpha_k x}{L}+\beta_k\bigg)\mathrm{d}x, v=\int_0^l \sin\bigg(\frac{\delta_m y}{l}+\theta_m\bigg)\sin\bigg(\frac{\delta_j y}{l}+\theta_j\bigg)\mathrm{d}y$$ For $$n=k, m=j$$, these integral evaluates to $$u=\frac{L}{2}-\frac{L}{4}[\sin(2\alpha_k + 2\beta_k)-\sin(2\beta_k)]$$, $$v=\frac{l}{2}-\frac{l}{4}[\sin(2\delta_k + 2\theta_k)-\sin(2\theta_k)]$$.

For $$n\neq k$$, the integrals are $$0$$ in this particular problem. I am skipping those details here.

Apart from these integrals, one also encounters the following while solving

$$I_1=\int_0^L \sin\bigg(\frac{\alpha_k x}{L}+\beta_k\bigg)=\frac{L}{\alpha_k}[\cos(\beta_k)-\cos(\alpha_k+\beta_k)], I_2=\int_0^l \sin\bigg(\frac{\delta_j y}{l}+\theta_j\bigg)=\frac{l}{\delta_j}[\cos(\theta_j)-\cos(\delta_j+\theta_j)]$$

I mention the $$u,v,I_1,I_2$$ expressions here as I guess these might help in simplification. I will also post their MMA code if someone wishes to use:

u = L/2 - (L/4)*(Sin[2 α + 2 β] - Sin[2 β])
v = l/2 - (l/4)*(Sin[2 δ + 2 θ] - Sin[2 θ])
I1 = (L/α)*(Cos[β] - Cos[α + β])
I2 = (l/δ)*(Cos[θ] - Cos[δ + θ])


Alternative approach This is an alternative approach from Andrea's great answer:

I tried an alternative approach where I write the two linear equations (i.e ortheq1 and ortheq2) containing $$C_1$$ and $$C_2$$ as the following:

$$A_{11}C_1+A_{12}C_2=XX_1 \tag4$$ $$A_{21}C_1+A_{22}C_2=XX_2 \tag5$$

I then extracted the coefficients of $$C_1$$ and $$C_2$$ using the following code, for which I took help from this answer given by Natas

Module[{coeffs = CoefficientList[Subtract @@ ortheq1, {C1, C2}]},
A11[C1_] = coeffs[[2, 1]] C1;
A12[C2_] = coeffs[[1, 2]] C2;
X11 = -coeffs[[1, 1]];
]
(Subtract @@ ortheq1) - (A11[C1] + A12[C2] - XX1) // Simplify
(* 0 *)

Module[{coeffs =
CoefficientList[Subtract @@ ortheq2, {C1, C2}]},
A21[C1_] = coeffs[[2, 1]] C1;
A22[C2_] = coeffs[[1, 2]] C2;
XX2 = -coeffs[[1, 1]];]

(Subtract @@ ortheq2) - (A21[C1] + A22[C2] - XX2) // Simplify
(*0*)


The solution of $$(4),(5)$$ is pretty straightforward in terms of $$A_{11},A_{12},A_{21},A_{22},XX_1$$ and $$XX_2$$

$$\text{C1}\to -\frac{A_{22} \text{XX}_1-A_{12} \text{XX}_2}{A_{12} A_{21}-A_{11} A_{22}},\text{C2}\to -\frac{A_{11} \text{XX}_2-A_{21} \text{XX}_1}{A_{12} A_{21}-A_{11} A_{22}}$$

I then simplified (using Mathematica and some hand calculations by looking at similar terms) $$A_{11},A_{12},A_{21},A_{22},XX_1$$ and $$XX_2$$ to get the following:

I ran your code and indeed, CC1 and CC2 are quite some expressions.

My approach to simplifying them is rather hacky, I cannot recommend doing this blindly - but here it goes:

1. Extract all subfactors, that only contain a big sum each, I have done this by hand for both expressions. For CC1:
C211 = CC1[][][];
C212 = CC1[][][];
C222 = CC1[][][];
C223 = CC1[][][];
C3111 = CC1[][][][];
C3112 = CC1[][][][];
C3122 = CC1[][][][];
C3123 = CC1[][][][];


and for CC2:

D12 = CC2[][];
D131 = CC2[][][];
D211 = CC2[][][];
D22 = CC2[][];
D2311 = CC2[][][][];
D2312 = CC2[][][][];
D2322 = CC2[][][][];
D2323 = CC2[][][][];
D24111 = CC2[][][][][];
D24112 = CC2[][][][][];
D24122 = CC2[][][][][];
D24123 = CC2[][][][][];

1. Let us take a look at the structure of the two expressions by replacing each of the summands with an appropriately named variable. For CC1 we get:
X1 = CC1 /. {C211 -> X211 , C212 -> X212, C222 -> X222, C223 -> X223,
C3111 -> X3111, C3112 -> X3112, C3122 -> X3122, C3123 -> X3123}


which gives $$-\frac{\text{X211} \text{X212}-\text{X222} \text{X223}}{\text{X212} \text{X3112}-\text{X223} \text{X3122}}$$,

and for CC2:

X2 = CC2 /. {D12 -> Q12, D313 -> Q313, D211 -> Q211, D22 -> Q22,
D2311 -> Q2311, D2312 -> Q2312, D2322 -> Q2322, D2323 -> Q2323,
D24111 -> Q24111, D24112 -> Q24112, D24122 -> Q24122,
D24123 -> Q24123}


which gives $$\frac{\text{Q22} (\text{Q211} \text{Q2311}-\text{Q12} \text{Q2323})}{\text{Q211} (\text{Q211} \text{Q24112}-\text{Q22} \text{Q2323})}-\frac{\text{Q12}}{\text{Q211}}$$.

This looks much friendlier. Also, notice that some of the $$Q$$-factors repeat, that means instead of 12 different factors we are actually only dealing with 6 different factors.

1. Simplify each factor while keeping in mind the overall expression. Obviously there is a lot of freedom of choice here. I strongly suggest you take a look at my suggestions and then implement something that makes more sense to you and the particular problem you are tackling.

I started with CC1 and simplified each of the six $$X$$-factors by themselves:

Y211 = Simplify[C211]
Y212 = Simplify[C212]
Y212 = FullSimplify[Y212]
Y222 = FullSimplify[C222]
Y223 = 1/(α (bc^2 + δ^2)^2)
FullSimplify[C223 /. {a_/(α (bc^2 + δ^2)^2) -> a}]
Y3111 = FullSimplify[C3111]
Y3112 = 1/(α (bc^2 + δ^2)^2)
FullSimplify[
Simplify[C3112 /. {a_/(α (bc^2 + δ^2)^2) -> a}]]
Y3122 = FullSimplify[C3122]
Y3123 = 1/(α (bc^2 + δ^2)^2)
FullSimplify[C3123 /. {a_/(α (bc^2 + δ^2)^2) -> a}]


Then I considered the structure of CC1 and calculated the nominator and denominator, simplifying both by substituting some of the factors of the denominator in the nominator. (essentially just $$\frac{a}{b} \to \frac{ca}{cb}$$, but I used the substitution since that is faster than multipliying and then simplifying):

R11 = (Y211 Y212 - Y222 Y223) /. { -1/(
8 α (bh^2 + α^2) δ^2 (bc^2 + \
δ^2)^2) -> (-64 α^2 (bh^2 + α^2)^2 δ^2 \
(bc^2 + δ^2)^2)/(
8 α (bh^2 + α^2) δ^2 (bc^2 + \
δ^2)^2) ,
1/(8 α^2 (bh^2 + α^2)^2 δ (bc^2 + \
δ^2)) -> (
64 α^2 (bh^2 + α^2)^2 δ^2 (bc^2 + \
δ^2)^2)/(
8 α^2 (bh^2 + α^2)^2 δ (bc^2 + \
δ^2))} /. E^(-bc - bh) l^2 L^2 ->  -1 /.
E^(-bc - bh + α_) l^2 L^2 ->  -E^α

R12 = (Y212 Y3112 - Y223 Y3122) /. -1/(
64 α^2 (bh^2 + α^2)^2 δ^2 (bc^2 + \
δ^2)^2) -> -1 /. -E^(-bc - bh + α_)
l^2 L^2 ->  -E^α


Now comes the most fidely bit - finding expressions that come up "often", giving them a new name and substituting. I did this one by one, overall I came up with:

    R112 = R11 /. -Sin[2 β] + Sin[2 (α + β)] -> uu /.
Sin[2 θ] - Sin[2 (δ + θ)] -> vv /.
Cos[β] - Cos[α + β] -> II1 /.
Cos[θ] - Cos[δ + θ] ->
II2 /. α Cos[α + β] +
bh Sin[α + β] ->
var1 /. -bh ph α Cos[α] + (bh^2 γ + \
α^2 (-ph + γ)) Sin[α] ->
var2 /. -bc pc δ Cos[δ] + (bc^2 γ + \
(-pc + γ) δ^2) Sin[δ] ->
var3 /. δ Cos[δ + θ] +
bc Sin[δ + θ] ->
var4 /. α Cos[β] - bh Sin[β] ->
var5 /. -α Cos[β] +
bh Sin[β] -> -var5 /. δ Cos[θ] -
bc Sin[θ] -> var6 /. -δ Cos[θ] +
bc Sin[θ] -> -var6 /.
2 α + Sin[2 β] - Sin[2 (α + β)] ->
var7 /. bh ph α Cos[α] + (bh^2 γ + \
α^2 (ph + γ)) Sin[α] -> var8  /.
bc pc δ Cos[δ] + (bc^2 γ + (pc + γ) \
δ^2) Sin[δ] -> var9

R113 = Simplify[R112] (*To see where I am at*)

(*I thought a second round might be good:*)
R113 /. bh ph var1 var5 α -> war1 /.
bc pc var4 var6 δ -> war2 /. -2 war1 +
E^bh (-ph α (bh^3 + (-1 +
bh) bh α^2 + α^4) + α (bh^2 + \
α^2)^2 γ -
var2 (bh^2 + α^2) Cos[α + 2 β]) ->
zar1 /.  δ (bc^3 (pc + bc ) +
bc ((-1 + bc) pc +
2 bc γ) δ^2 + (pc + γ) δ^4) \
-> zar2 /. bc^3 (-pc + bc γ) δ +
bc (pc - bc pc +
2 bc γ) δ^3 + (-pc + γ) δ^5 ->
zar3 /. bc^2 + δ^2 -> rad1^2 /.
bh^2 + α^2 ->
rad2^2 /. δ (bc^3 (pc + bc γ) +
bc ((-1 + bc) pc +
2 bc γ) δ^2 + (pc + γ) δ^4) ->
zar4 /. ph α (bh^3 + (-1 + bh) bh α^2 + α^4) +


and the same goes for R12:

    R122 = R12 /. -Sin[2 β] + Sin[2 (α + β)] -> uu /.
Sin[2 θ] - Sin[2 (δ + θ)] -> vv /.
Cos[β] - Cos[α + β] -> II1 /.
Cos[θ] - Cos[δ + θ] ->
II2 /. α Cos[α + β] +
bh Sin[α + β] ->
var1 /. -bh ph α Cos[α] + (bh^2 γ + \
α^2 (-ph + γ)) Sin[α] ->
var2 /. -bc pc δ Cos[δ] + (bc^2 γ + \
(-pc + γ) δ^2) Sin[δ] ->
var3 /. δ Cos[δ + θ] +
bc Sin[δ + θ] ->
var4 /. α Cos[β] - bh Sin[β] ->
var5 /. -α Cos[β] +
bh Sin[β] -> -var5 /. δ Cos[θ] -
bc Sin[θ] -> var6 /. -δ Cos[θ] +
bc Sin[θ] -> -var6 /.
2 α + Sin[2 β] - Sin[2 (α + β)] ->
var7 /. bh ph α Cos[α] + (bh^2 γ + \
α^2 (ph + γ)) Sin[α] -> var8  /.
bc pc δ Cos[δ] + (bc^2 γ + (pc + γ) \
δ^2) Sin[δ] -> var9

R123 = Simplify[R122]

R123 /. bh ph var1 var5 α -> war1 /.
bc pc var4 var6 δ -> war2 /. -2 war1 +
E^bh (-ph α (bh^3 + (-1 +
bh) bh α^2 + α^4) + α (bh^2 + \
α^2)^2 γ -
var2 (bh^2 + α^2) Cos[α + 2 β]) ->
zar1 /.  δ (bc^3 (pc + bc γ) +
bc ((-1 + bc) pc +
2 bc γ) δ^2 + (pc + γ) δ^4) \
-> zar2 /. bc^3 (-pc + bc γ) δ +
bc (pc - bc pc +
2 bc γ) δ^3 + (-pc + γ) δ^5 ->
zar3 /. bc^2 + δ^2 -> rad1^2 /.
bh^2 + α^2 ->
rad2^2 /. δ (bc^3 (pc + bc γ) +
bc ((-1 + bc) pc +
2 bc γ) δ^2 + (pc + γ) δ^4) ->
zar4 /. ph α (bh^3 + (-1 + bh) bh α^2 + α^4) +


Overall, this allowed me to find $$CC1 = \frac{A1-A2}{A3}$$:

    (*"Final" Result for CC1 = R1*)
A1 = 8 II2 ph rad2^2 (Ta - thi) (uu - 2 α) α (2 war2 +
E^bc zar4 -
2 θ]) (-α Cos[α + β] +
E^bh (α Cos[β] + bh Sin[β]) -
bh Sin[α + β]);
A2 = 8 E^(-w γ)
II1 pc rad1^2 (Ta - tci) zar1 δ (vv +
2 δ) (E^
bc δ Cos[θ] - δ Cos[δ + θ] +
bc E^bc Sin[θ] - bc Sin[δ + θ]);
A3 = E^(-w γ) (vv +
2 δ) (-E^(
2 w γ) (uu - 2 α) (2 war1 +
E^bh (zar5 - rad2^2 var8 Cos[α + 2 β])) (2 war2 +
E^bc zar2 - E^bc rad1^2 var9 Cos[δ + 2 θ]) -
var7 zar1 (-2 war2 +
E^bc (zar3 - rad1^2 var3 Cos[δ + 2 θ])));
R1 = ( A1 - A2)/A3


Still not pretty, but certainly better than what we started with.

For CC2 I followed the same steps and of course tried to reuse the substitutions I already made for CC1.

• This is some awesome work Andrea. Thank you and much appreciation. This must have taken quite some time and effort. From the time I asked this question, I had tried an alternative approach which I have added to my original question, have a look if you are interested. Nevertheless, you have really managed to get the constants quite simplified. Jul 8 '20 at 17:05