4
$\begingroup$

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[[1]] - bc1[[2]])*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[[1]] - bc2[[2]])*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:

enter image description here enter image description here

$\endgroup$
2
+50
$\begingroup$

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[[2]][[1]][[1]];
C212 = CC1[[2]][[1]][[2]];
C222 = CC1[[2]][[2]][[2]];
C223 = CC1[[2]][[2]][[3]];
C3111 = CC1[[3]][[1]][[1]][[1]];
C3112 = CC1[[3]][[1]][[1]][[2]];
C3122 = CC1[[3]][[1]][[2]][[2]];
C3123 = CC1[[3]][[1]][[2]][[3]];

and for CC2:

D12 = CC2[[1]][[2]];
D131 = CC2[[1]][[3]][[1]];
D211 = CC2[[2]][[1]][[1]];
D22 = CC2[[2]][[2]];
D2311 = CC2[[2]][[3]][[1]][[1]];
D2312 = CC2[[2]][[3]][[1]][[2]];
D2322 = CC2[[2]][[3]][[2]][[2]];
D2323 = CC2[[2]][[3]][[2]][[3]];
D24111 = CC2[[2]][[4]][[1]][[1]][[1]];
D24112 = CC2[[2]][[4]][[1]][[1]][[2]];
D24122 = CC2[[2]][[4]][[1]][[2]][[2]];
D24123 = CC2[[2]][[4]][[1]][[2]][[3]];
  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.

The interesting task is now

  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) +
        rad2^4 α γ -> zar5

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) +
        rad2^4 α γ -> zar5

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 - 
         E^bc rad1^2 var9 Cos[δ + 
            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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Indrasis Mitra Jul 8 at 17:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.