User:Egm6936.f10/Stochastic 1-D heat problem

In this section, we will employ a stochastic 1-D heat conduction example to illustrate some concepts used in gPC to handle stochastic problems. This example mainly comes from D. Xiu et al. 2002. The formulation is basically developed from Prof. Loc Vu-Quoc's lecture notes [[media:pea2.s12.sec60.a.djvu| Principle of Engineering Analysis (Spring2012)_Sec.60.a]]. For any other information, one should resort to the original articles.

=Stochastic 1-D heat problem= The statement of the problem is as follows: The boundary conditions are and the random diffusivity is where $$\displaystyle \omega$$ is the event, outcome, $$\displaystyle \epsilon(\omega)$$ is the associated random variable, and $$\displaystyle \kappa(x;\omega)>0$$. Before starting to solve this problem, we provide the exact solution of this problem for comparison with the gPC solution. The exact solution is For $$\displaystyle \epsilon(\omega)\neq 0:$$ For $$\displaystyle \epsilon(\omega)= 0:$$

=Legendre Chaos and Uniform distribution=

In this subsection, we assume 

(see Section Note for the meaning of $$\sigma$$), where

is a continuous uniform random variable. In D. Xiu et al. 2002, there is not much information on $$\displaystyle \sigma $$. But in another paper, D. Xiu et al. 2003, Xiu termed $$\displaystyle \sigma $$ as the standard deviation. In fact, $$\displaystyle \sigma $$ is a scaling factor that is used to scale the standard deviation of the uniform variable, since for a given Uniform distribution, the standard deviation is known; see Uniform distribution. It also can be verified using the definition of standard deviation(see Section Note for the meaning of $$\sigma$$). The polynomial chaos associated with uniform distribution is the Legendre polynomial chaos.

From stochastic to deterministic
Applying Legendre Chaos expansion to $$\displaystyle \kappa (x;\omega)$$ and $$\displaystyle u (x;\omega)$$, we have

Generally, the expansion is infinite. But in practice, a finite expansion is adopted, and when to truncate the expansion is problem dependent. Here, the summations are truncated at $$\displaystyle N_{\kappa}$$ and $$\displaystyle N_u$$ respectively. And it is not necessary that all expansions should have the same truncation number.

Truncation number $$\displaystyle N_{\kappa}$$ and $$\displaystyle N_u$$ also are the highest order of the employed polynomials in each expansion. The random parameter $$\displaystyle \omega$$ is absorbed into the Legendre basis $$\displaystyle P(\xi(\omega))$$. For simplicity, the notation $$\displaystyle \xi$$ and $$\displaystyle P(\xi)$$ will be used in following content unless confusion arise. The coefficients $$\displaystyle \kappa_{i} $$ and $$\displaystyle u_{j} $$ are deterministic. Plugging the Legendre expansion (ab6b) into governing equation (ab1), we obtain Upon simplification, equation (ab7) can be rewritten as A Galerkin projection of the above equation onto each Legendre polynomial basis $$\displaystyle \{P_{k}\}$$ is then conducted to ensure that the error is orthogonal to the Legendre polynomial space. By projecting and employing the orthogonality property of Legendre basis, we obtain for each $$\displaystyle k=0,1,...,N_{u}$$.

Or written in the matrix form as: 

In above equation (ab9a), $$\displaystyle c_{ijk} = \langle P_{i}P_{j}P_{k}\rangle$$, which is the weighted inner product, can be independently evaluated from the Legendre polynomial basis. The relation between number of final equations $$\displaystyle P+1 $$, the dimensionality of the random inputs $$\displaystyle M $$ and the highest order of the polynomials $$\displaystyle N_{u}$$ is $$\displaystyle (P+1) = \frac{(M+N_u)!}{(M)!(N_u)!} $$. In this case, the dimensionality of the random inputs $$\displaystyle M = 1 $$. Hence, equation (ab9b) is a set of $$\displaystyle (N_u+1)$$ coupled PDEs.

Zeros of C matrix
Before proceeding to derivation of the governing equation system, it's necessary to look at the zero entries of the coefficient matrix

that appears in Eq.(ab9b).

Oddness and evenness of Legendre polynomial
For $$\displaystyle k=0,1,2,...$$, Legendre polynomial $$\displaystyle P_{2k}$$ is even and $$\displaystyle P_{2k+1}$$ is odd.

Above statement can be easily verified using the explicit form of Legendre polynomial:

For more details about oddness and evenness of Legendre polynomial, please see lecture notes [[media:pea2.s12.sec43.a.djvu| Principle of Engineering Analysis (Spring2012)_Sec.43.a]]. Also, some basic properties of odd and even function are listed below:
 * The product of two even functions is even.
 * The product of two odd functions is even.
 * The product of an even function and an odd function is odd.
 * The integral of an odd function over skew-symmetric domain is zero.

Zero coefficients in C matrix
As by the definition, the coefficient matrix $$\displaystyle \mathbf C_k $$ is a integral of triple product of Legendre polynomials over $$\displaystyle [-1,1] $$.

Using those basic properties of odd and even functions (i.e., the function parity), we can arrive at the following 2 criteria for deciding the zero entries of the $$\displaystyle \mathbf C_k $$ matrix:

A.1) For coefficient $$\displaystyle c_{ijk} $$ with nonzero indices,  if the sum $$\displaystyle (i+j+k)$$ is an odd number, then this $$\displaystyle c_{ijk} = 0 $$.

For the sum $$\displaystyle (i+j+k)$$ to be an odd number,

we have 2 cases: (a) Only one of the three indices is odd (and the other two indices are even), and (b) all three indices are odd.

For both the above two cases, the resulting integrand is an odd function, and thus the integration of this odd function over a symmetric domain is zero.

 A.2) For the coefficients $$\displaystyle c_{ijk} $$ in which one index is zero, e.g., say $$\displaystyle k=0$$ and $$\displaystyle P_k (x) = P_0 (x) = 1$$, then we only need to look at the product $$\displaystyle P_i P_j$$ and use the orthogonality of the Legendre polynomials; there are 2 cases: (a) if $$\displaystyle i=j$$, then $$\displaystyle c_{ii0} \ne 0$$; (b) if $$\displaystyle i \ne j$$, then $$\displaystyle c_{ij0} = 0$$.

For example, consider $$\displaystyle c_{240}$$. Although $$\displaystyle 2+4+0 = 6$$ is even, we have $$\displaystyle c_{240} = 0$$ due to the second condition.

Gupta and Narasimhan (2007) wrote the above 2 criteria in more mathematical form as follows:

The coefficient $$\displaystyle c_{ijk} = 0 $$, if either of following two conditions is satisfied:

B.1) the remainder $$\displaystyle (i + j + k) \text{ mod } 2 $$ is nonzero.

B.2) the combination of indices $$\displaystyle i, j, k $$ doesn't satisfy the triangle inequality, i.e., $$\displaystyle i + j < k $$ or $$\displaystyle j + k < i $$ or $$\displaystyle i + k < j $$.

The above 2 conditions can be explained as follows.

Note that $$\displaystyle (i + j + k) \text{ mod } 2 $$ means the remainder of the division of $$\displaystyle (i + j + k)$$ by 2.

The first condition can be rewritten as

\begin{cases} (i+j+k) \text{ mod } 2 = 0 & \text{then } c_{ijk} \ne 0 \\ (i+j+k) \text{ mod } 2 = 1 & \text{then } c_{ijk} = 0 \end{cases} $$

For example, consider $$\displaystyle c_{140}$$; since $$\displaystyle 1+4+0 = 5$$ is odd, we have $$\displaystyle c_{140} = 0$$.

Another example, since $$\displaystyle 2+4+0 = 6$$ is even, we have $$\displaystyle c_{240} = 0$$.

The reason is that since one index is zero, say $$\displaystyle k=0$$, we need to see whether the other 2 indices $$\displaystyle i,j$$ are equal or not.

If $$\displaystyle i=j$$, then $$\displaystyle c_{ii0} \ne 0$$;

If $$\displaystyle i \ne j$$, then $$\displaystyle c_{ij0} = 0$$ due to the orthogonality of the Legendre polynomials, as already explained above.

For the 2nd condition, i.e., the triangle inequalities, it's more general than the case in which one index is zero.

For example, $$\displaystyle c_{237} = 0$$, we can't conclude this result from any of the first two conditions, A.1 and A.2, but from the triangle inequalities, B.2.

MATLAB symbolic integration function int is used to check the result.



The above two conditions can be used to spot all zero coefficients in the $$\displaystyle \mathbf C_k $$ matrix.

Since the indices in matrix coefficients in matlab begin with 1 instead of with 0, we introduce the following notation $$\displaystyle \bar {\mathbf C}_r := [ \bar c_{pqr} ] $$ with $$\displaystyle p,q,r = 1,...,(N+1)$$ so that $$\displaystyle \bar c_{pqr} = c_{(p-1)(q-1)(r-1)}$$.

Below is a figure showing the matrix profiles (in which non-zero coefficients are represented as blue dots) of the matrices $$\displaystyle \bar {\mathbf C}_r := [ \bar c_{pqr} ]$$, for $$\displaystyle r=1,...,(N+1)$$, with $$\displaystyle N = N_{\kappa} = N_{u} = 7$$.

The matrix profiles are obtained with the MATLAB function spy.

In the matrix-profile figure, the following legend is used:

The case r=1 or k=0 corresponds to the matrix $$\displaystyle \bar {\mathbf C}_1 = \mathbf C_0$$ in Eq.(ab57).

The case r=2 or k=1 corresponds to the matrix $$\displaystyle \bar {\mathbf C}_2 = \mathbf C_1$$ in Eq.(ab58), where it can be seen that the diagonal contains only zeros, and the two off-diagonals have non-zero coefficients (these 2 off-diagonals correspond to the 2 off-diagonals with blue dots in the matrix profile).

and so on so forth.



In order to demonstrate the efficiency of the two conditions, B.1 and B.2, made on the calculation of the triple product of Legendre polynomials, we provide one more example with a larger number of Legendre polynomials, $$\displaystyle N = 19$$. The MATLAB functions tic and toc are used to calculate the computation time for both codes, i.e, the one with the conditions and the other without.

For the code with conditions B.1 and B.2, the elapsed time is 180.983425 seconds (for computing $$\displaystyle \mathbf C_{0}, \cdots , \mathbf C_{19}$$ ),

while for the one without these conditions, the elapsed time is 870.205230 seconds.

As can be seen in above computation time, the two conditions greatly improve the computing efficiency. As the matrix size becomes even more larger, say over one hundred, these two conditions are more and more important in saving computing time.

Three matrix profiles for $$\displaystyle \bar {\mathbf C}_{1} = \mathbf C_{0}$$, $$\displaystyle \bar {\mathbf C}_{10} = \mathbf C_{9}$$ and $$\displaystyle \bar {\mathbf C}_{20} = \mathbf C_{19}$$, and for $$\displaystyle N = 19$$, are provided below:



Governing equation systems
In this section, we will elaborate on how to determinize the governing equation for both non-uniform truncation and uniform truncation cases. First, let's look at the uniform case.

Uniform truncation
For simplicity, we will use the same truncation number $$\displaystyle N$$ for both expansions of $$\displaystyle \kappa$$ and $$\displaystyle u$$ in following content. We will write out the equations for different $$\displaystyle N (N = 1, 2, ..., 7)$$ to show how these equations will look like. And later, some details on how to the use spectral method to solve the equation system is given.

N = 1
When $$\displaystyle N_{\kappa} = N_{u} = N = 1 :$$ The Legendre Chaos bases are Therefore, we can immediately compute the coefficients $$\displaystyle c$$ as follows:     And the two equations are:   Since the conductivity is already given as $$\displaystyle \kappa(x;\omega)=1+\sigma \xi(\omega)x $$, we can directly have the components of $$\displaystyle \kappa_{i}(x), i = 0,1$$ as Substituting above quantities into equation (ab15) and (ab16), these two equations can be further modified as

N = 2
When $$\displaystyle N_{\kappa} = N_{u} = N = 2 :$$ The Legendre Chaos bases are Except above calculated values for the coefficient $$\displaystyle c$$, we have additionals as follows:

In this case, we have $$\displaystyle \kappa_{2}(x) = 0$$.

Substituting above known quantities into equation (ab9), the three equations can be simplified as follows:

N = 3
When $$\displaystyle N_{\kappa} = N_{u} = N = 3 :$$

The Legendre Chaos bases are In this case, we have $$\displaystyle \kappa_{3}(x) = 0$$. The details of the coefficients $$\displaystyle c$$ will not given in this and following subsections, instead, a MATLAB code for calculating these quantities will be provided in the end of this section. Finally, the four equations can be simplified as follows:

N = 4
When $$\displaystyle N_{\kappa} = N_{u} = N = 4 :$$

The Legendre Chaos bases are

In this case, we have $$\displaystyle \kappa_{4}(x) = 0$$.

Finally, the five equations can be simplified as follows:

N = 5
When $$\displaystyle N_{\kappa} = N_{u} = N = 5 :$$

The Legendre Chaos bases are

In this case, we have $$\displaystyle \kappa_{5}(x) = 0$$.

Finally, the six equations can be simplified as follows:

N = 6
When $$\displaystyle N_{\kappa} = N_{u} = N = 6 :$$ The Legendre Chaos bases are

In this case, we have $$\displaystyle \kappa_{6}(x) = 0$$.

Finally, the six equations can be simplified as follows:

N = 7
When $$\displaystyle N_{\kappa} = N_{u} = N = 7 :$$ The Legendre Chaos bases are

In this case, we have $$\displaystyle \kappa_{7}(x) = 0$$.

Finally, the six equations can be simplified as follows:

These equations are deterministic and can be solved by any conventional method, e.g., finite difference method, finite element method, etc.

Non-uniform truncation
In this section, we use different truncation number for $$\displaystyle \kappa$$ and $$\displaystyle u$$. Since the focus is to find how the result will different from the uniform truncation case, we will only provide the case when $$\displaystyle N_{\kappa} = 3$$ and $$\displaystyle N_{u} = 4$$ here.  $$\displaystyle N_{\kappa} = 3, N_{u} = 4 $$ : The Legendre Chaos bases are The coefficients $$\displaystyle \kappa_{i}(x), i = 0,1,2,3$$ are

The five governing equations can be simplified as follows:

Why the last column entries of above $$\displaystyle \mathbf C$$ are all zeros?

The $$\displaystyle \mathbf C$$ matrix for this equation is:

According to the Condition 2 in the section of Zero entries of the $$\displaystyle \mathbf C$$, it's obviously the last column entries should all be zeros.

The boundary condition for this case is the same as that when uniform $$\displaystyle N_{\kappa}=N_u=N = 4$$ is used. See Boundary condition section N=4 for more details.

Note: As you can observe from above five deterministic governing equations, due to the property that the conductivity $$\kappa$$ can be exactly represented by the first two Legendre polynomials, the final form of these equations won't be changed given more series expansion for the temperature $$u$$. Rather than in this 1-D heat conduction problem, the solution of non-uniform truncation should be different from uniform truncation.

Boundary conditions for uniform truncation case
The way to determinize the boundary conditions is the same as that used for the governing equation.

N = 1
When $$\displaystyle N_{\kappa}=N_u=N = 1 :$$

The stochastic Galerkin's method is employed here to get rid of the uncertainties in equation(ab65).

where $$\displaystyle a_{ij}= \mathbb E \langle P_i P_j \rangle$$ is the weighted inner product, and $$\displaystyle p_i = \mathbb E \langle P_i 1 \rangle = \mathbb E \langle P_i P_0 \rangle $$. The value of $$\displaystyle a_{ij}$$ and $$\displaystyle p_i$$ can be evaluated using the orthogonal property of Legendre polynomials.

Plugging above coefficients into equation(ab66), the boundary conditions are:

N = 2
When $$\displaystyle N_{\kappa}=N_u=N = 2 :$$

The same procedure is used for $$\displaystyle N = 2 $$.

Finally, we have the boundary conditions as

N = 3
When $$\displaystyle N_{\kappa}=N_u=N = 3 :$$ The essential boundary conditions are as following

N = 4
 When $$\displaystyle N_{\kappa}=N_u=N = 4 :$$ The essential boundary conditions are as following

N = 5
When $$\displaystyle N_{\kappa}=N_u=N = 5 :$$ The essential boundary conditions are as following

N = 6
When $$\displaystyle N_{\kappa}=N_u=N = 6 :$$ The essential boundary conditions are as following

N = 7
When $$\displaystyle N_{\kappa}=N_u=N = 7 :$$ The essential boundary conditions are as following

Solution of deterministics
For simplicity, only the solution When N = 2 is provided.

$$\displaystyle \sigma = 0.1 $$

$$\displaystyle \sigma = 0.5 $$

$$\displaystyle \sigma = 0.9 $$

The mean-square error, $$\displaystyle e_2(x)$$, of the numerical solution from the gPC expansion $$\displaystyle u_p(x,\omega)$$ is computed $$\displaystyle e_2(x) = (\mathbb E[u_p(x,\omega)-u_e(x,\omega)]^2)^{(1/2)}$$ The $$\displaystyle L_\infty$$ norm of $$\displaystyle e_2(x)$$ is plotted with respect to $$\displaystyle p$$, the order of polynomial chaos expansion.



=Note= $$\displaystyle 1\bullet$$ The meaning of $$\displaystyle \sigma $$ in (Eq.(ab5b)): As mentioned in previous section, that $$\displaystyle \sigma $$ is a scaling factor applied to the standard deviation of the random variable $$\displaystyle \xi $$ which assumed to have a uniform distribution. According to D. Xiu et al. 2002 and D. Xiu et al. 2003, $$\displaystyle \sigma $$ is used to evaluate the performance of Polynomial Chaos Expansion at different 'distribution density', and also compare the results with perturbation method, which is valid only when $$\displaystyle \sigma $$ is very small, say $$\displaystyle \sigma = 0.1$$. By definition, the standard deviation of a random variable $$\displaystyle X $$, is the square root of its variance. The variance of $$\displaystyle X $$ is given by

where $$\displaystyle \mu $$ is the expected value (mean) of $$\displaystyle X $$. For $$\displaystyle X \doteq \mathcal{U}\left ( a,b \right )$$, the mean of random variable $$\displaystyle X $$ is Hence, the standard deviation of random variable $$\displaystyle X $$ is Notice that, generally, $$\displaystyle \sigma $$ is adopted to notate the standard deviation, while at the beginning, we have used $$\displaystyle \sigma $$ to denote a scaling factor applied to the random variable $$\displaystyle \xi $$. According to above formula, the standard deviation for $\displaystyle \epsilon $, ($$\displaystyle \epsilon(\omega) = \sigma \xi(\omega)$$, where $$\displaystyle \xi (\omega)\doteq \mathcal{U}\left ( -1,1 \right )$$ is a continuous uniform random variable), is

$$\displaystyle 2\bullet$$ The components of $$\displaystyle \mathbf C_k $$ Matrix: In Eq.(ab9b), Galerkin's weighted residual form is used to handle the uncertain part in the polynomial chaos expansion. While doing projection, the weight can be both density associated with the polynomial chaos basis, here is $$\displaystyle \frac{1}{2} $$ for Legendre Polynomials, and unity, which will not results in any inconsistency since the equation still being balanced by multiplying any constant, except zero, on both sides. In above derivation, the density function $$\displaystyle \frac{1}{2} $$ is used, as you can see from Eq.(ab11). Let's check whether it will affect the solution of the equation system for the case when $$\displaystyle N =1 $$. Not including the associated density function for the Legendre polynomials in our case, the components of matrix $$\displaystyle \mathbf C_k $$ can be calculated as Comparing above computation of the $$\displaystyle \mathbf C_k $$ components with Eq.(ab11), Eq.(ab12), Eq.(ab13) and Eq.(ab14), we can see the difference between using and not using density function in the calculation of $$\displaystyle \mathbf C_k $$ matrix. The resulting equations are: Above two equations absolutely are the same as Eq.(ab15) and Eq.(ab16).

$$\displaystyle 3\bullet$$ Some coding issues: Symbolic operation was first adopted when the code to solve the coupled equation system was developed. This handles well for cases when $$\displaystyle \sigma = 0.5 $$ and $$\displaystyle \sigma =0.9 $$, and for only two, or less, polynomial chaos bases are employed in the expansion for case of $$\displaystyle \sigma = 0.1 $$. But when more polynomial chaos bases were used, numerical issues, that double function can't convert the symbolic value to numerical, rendered the final $$\displaystyle \mathbf L_2 $$ evaluation becomes unavailable. By means of using numerical operation to calculate the $$\displaystyle \mathbf L_2 $$ and $$\displaystyle \mathbf L_\infty $$ norms, this issue is well solved.

After running above code, following error message will pop out: ??? The following error occurred converting from sym to double: Error using ==> mupadmex Error in MuPAD command: DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use the VPA function instead. This error message will not disappear even after using VPA function instead. As far as I can think of, the reason should be the double function can't convert the symbolic expression into numerical number, since the expression is too large which beyond the limit of the double function.

Following is the symbolic expression:

=References=

=References=