According to Wikipedia, computer vision is an interdisciplinary field that deals with how computers can be made for gaining high-level understanding from digital images or videos. Computer vision tasks include methods for acquiring, processing, analyzing and understanding digital images, and extraction of high-dimensional data from the real world in order to produce numerical or symbolic information. It involves the development of a theoretical and algorithmic basis to achieve automatic visual understanding. As a scientific discipline, computer vision is concerned with the theory behind artificial systems that extract information from images. The image data can take many forms, such as video sequences, views from multiple cameras. In my graduate studies at Carnegie Mellon University, I took a course in computer vision, which is my favorite course of all time as well as the most difficult! These are projects I completed that will encapsulate the basic idea of this powerful computer science concept. Please note that I cannot provide my source code because these assignments may be reused in the course.

Tracking

The first groundbreaking work on template tracking was the Lucas-Kanade tracker. It basically assumes that the template undergoes constant motion in a small region. The Lucas-Kanade Tracker works on two frames at a time, and does not assume any statistical motion model throughout the sequence. The algorithm estimates the deformations between two image frames under the assumption that the intensity of the objects has not changed significantly betweenthe two frames. Starting with a rectangle $R_t$ on frame $I_t$, the Lucas-Kanade Tracker aims to move it by an offset (u; v) to obtain another rectangle $R_{t+1}$ on frame $I_{t+1}$, so that the pixel squared difference in the two rectangles is minimized. According to Wikipedia, $A^TA$ is the structure tensor of the image at all points in the rectangle. It's a gradient matrix that summarizes the principal directions in a specified neighborhood of that point.

$$A=\left[ \begin{array}{cc} I_x(q_1) & I_x(q_1) \\ I_x(q_2) & I_x(q_2) \\ \vdots & \vdots \\ I_x(q_n) & I_x(q_n) \\ \end{array} \right] \rightarrow A^TA = \left[ \begin{array}{cc} \sum I_x^2 & \sum I_xI_y \\ \sum I_xI_y & \sum I_y^2 \\ \end{array} \right] $$

As for the conditions, $A^TA$ should be invertible, well-conditioned and should not be too small due to noise (meaning the eigenvalues should not be too small). Below is a tracking example of a car with little variations. This is totally awesome because I always see bank robberies and car pursuit on T.V sky cams and now I got to learn the technology and science behind it!

The above video is quite simple but real data is often corrupted by unknown image noise or under varying illumination conditions. One way to address issue of appearance variation is to use the principal component analysis.

As you can see, this example isn't so perfect because of the variation in appearance such as shading and lighting. Can you think of a way to make it "easier" for computers to track these objects? How about reducing the color channels by converting the frames to grayscale? Or use various empirical affine motion equations to form a linear combination of bases.

Spatial Pyramid Matching for Scene Classification

The question of interest for this topic is "Given an image, can a computer program determine where it was taken?" I used the Bag-of-Words (BoW) approach. The BoW essentially trains the computer to reprense the world with visual words. The idea is to give the computer many(!) data (images) to understand and recognize so that it is "smart" enough to identify the next datum. This is where data science like Machine Learning and Data Mining come to play to supervise the dataset. To help the computer interpret the data, we filter the images and extract their properties. A filter bank is a series of flters that captures different visual properties. The Gaussian blurs are to smooth out the details by reducing noise. For example, convolving with the Laplacian Gaussian kernel allows easier detection of edges. Furthermore, the derivative of the Gaussian on the x-direction detects a change in intensity along that axis, thus vertical edges are more prominent. The same can be said for gradient of the Gaussian in respect to the y-axis which shows details on the horizontal edges. Bag of words is simple and efficient, but it discards information about the spatial structure of the image and this information is often valuable. One way to alleviate this issue is to use spatial pyramid matching. The general idea is to divide the image into a small number of cells, and concatenate the histogram of each of these cells to the histogram of the original image, with a suitable weight. The following outlines the basic idea of the algorithm.

1. Create a filter bank with different sort of scales, orientations, and bandwidths to measure the responses of the images.
2. Randomly select a number of pixels to examine and group these pixels' response in a cluster. From these data we are able to compile a dictionary of visual words.
3. Run through all images and obtain the visual word with the minimum distance on each pixel. Extract the frequency histogram where the number of bins is the number of words. For a more accurate result, use Spatial Pyramid Matching
4. Find the instance word with the largest similarity and assign its category to the image.

Keypoints - Detectors, Descriptors and Matching

Interest point detectors find particularly salient points in an image upon which we can extract a feature descriptor. In our case, we will be using BRIEF. Once we have extracted the interest points, we can use descriptors to match them between images to do neat things like panorama stitching or scene reconstruction. Keypoints are found by using the Difference of Gaussian (DoG) detector. This detector finds points that are extrema in both scale and space of a DoG pyramid. The objective of this assignment was to learn a tomato soup can and be able to recognize and match it from various angles. We accomplish this by identifying points of interest and match them with a descriptor. One simple descriptor is BRIEF which uses the Hamming distance to compute the similarity between the potential matches. The belowed image on the left shows the ideals points of critical contrast that I needed to match and my results on the right image.

Unfortunately, my algorithm does not yield any correct matches. The closest and best result I obtained was the top left. Of course, if the test image compares to itself, the percentage of correct matches would be 100%, however I believe once rotated, the derivative of Gaussian is more difficult in terms of computing the gradient at the edges. In other words, the eigenvalues is distorted due to the change of orientation.

Homographies

Robots often deal with planes, whether in the form of walls, ground, or some other at surface. When two cameras observe a plane, there exists a relationship between the images captured. This relationship is defined by a 3x3 transformation matrix, called a planar homography. This was definitely one of my more favorite topics to learn in Computer Vision because it is highly practical. Let's suppose a robot is at war, it needs to be able to "stitch" the planar images that it "sees" to interpret it as a warzone and take action. A planar homopgrahy allows us to compute how a planar scene would look from a second camera location give only the first image! Furthermore, we can extrapolate any camera angle from any location without knowing any internal camera parameters. Can you quess what simple use of homography? Panoramas! For this assignment, I was given two images of the Taj Mahal and I had to use homographies to create a panorama image of the same scene.

I designed and built a simple software package/interface for finite element analysis. The main idea of this project was to provide students access to analytical studies of the finite element method. Students should be able to obtain a basic understanding of the importance of FEM and how it works at the elementary level. Note that some calculus (including linear algebra and differential equations) is required to understand the mathematics.

Click on the icon to download the graphical user interface.

Click on the icon to download the lesson PDF.

Linear & Quadratic Interpolation of Continuous Function

§ Methods of Weighted Residual to Obtain Governing Weak Form

Consider: $-(k\hat{u}')'+b\hat{u}'+c\hat{u}=f$
Residual: $r(\hat{u})=-(k\hat{u}')'+b\hat{u}'+c\hat{u}-f$
Multiply by test function v such that v and u are in the same space (thus v also satisfies EBC) and integrate over domain: $$\int_0^L r(\hat{u})vdx = 0 \rightarrow \int_0^L [-(k\hat{u}')'+b\hat{u}'+c\hat{u}-f]vdx=0 $$ Suppose v is sufficiently smooth so that we can integrate the 1st term by parts: $$\int_0^L -(k\hat{u}')'vdx=-(k\hat{u}')v|^L_0 + \int_0^Lk\hat{u}'v'dx$$ Suppose the initial boundary is zero, then u(0)=0, it follows that that v at x=0 must also be 0 and we have $$\int_0^L k\hat{u}'v'dx + \int_0^Lc\hat{u}vdx - \int_0^Lfvdx - T_vv(L)=0$$

§ Numerical Approximation: Galerkin Method in Weak Form

Let $\hat{u} \sim u_N = \sum_{k=1}^N \alpha_k \phi_k(x)$, then $\hat{u} \sim u_N = \sum_{i=1}^N \beta_i \phi_i(x)$ $$\int_0^L k \left( \sum_{k=1}^N \alpha_k \phi_k'(x) \right)\left( \sum_{i=1}^N \beta_i \phi_i'(x) \right)dx + \int_0^L c \left( \sum_{k=1}^N \alpha_k \phi_k(x) \right)\left( \sum_{i=1}^N \beta_i \phi_i(x) \right)dx - \int_0^L f \left( \sum_{i=1}^N \beta_i \phi_i(x) \right)dx-T_v \left( \sum_{i=1}^N \beta_i \phi_i(L) \right)dx$$ Let $$K_{kj}=\int_0^L k(x)\phi_j'(x)\phi_k'(x)dx=K_{jk}$$ $$C_{kj}=\int_0^Lc(x)\phi_j(x)\phi_k(x)dx$$ $$f_k=\int_0^L f\phi_k(x)dx+\phi_k(L)$$ The final equation to solve for is $$\sum_i B_i \left\{ \sum_k (K_{jk}+C_{jk})a_k -f_i \right\}=0 \rightarrow \left\{ \sum_k (K_{jk}+C_{jk})a_k -f_i \right\}=0 \forall B_i, \quad i=1, \dotsc ,n$$ The goal is to solve for $\alpha_k$ in $D\alpha=f$ where $D_{jk}=K_{jk}+C_{jk}$

§ Selecting Appropriate Basic Functions

Consider the following linear basis functions

Trial functions: $U_n(x)=\sum_{i=1}^n \alpha_i \phi_i(x)$
Test functions: $V_n(x)=\sum_{i=1}^n \beta_i X_i(x)$
How do we choose a $\phi(x)$ such that $\hat{u}$ satisfies the boundary conditions?

If we consider using quadratic basis

Local Functions: $$\left\{ \begin{array}{c} \phi_A^e(\xi)=-\frac{1}{2}(\xi)(1-\xi)\\ \phi_B^e(\xi)=\frac{1}{2}(\xi)(1+\xi)\\ \phi_C^e(\xi)=(1-xi)(1+\xi)\\ \end{array} \right.$$ If we notice that $x=x_c+\xi \frac{h}{2} \rightarrow \xi=\frac{2(x-x_c)}{2}$ $$\left\{ \begin{array}{c} \phi_A^e(\xi)=\frac{-x(x_C)}{h} \frac{h-2x+2x_C)}{h}\\ \phi_B^e(\xi)=\frac{x(x_C)}{h} \frac{h+2x-2x_C)}{h}\\ \phi_C^e(\xi)=\frac{h-2x+2x_C}{h} \frac{h+2x+2x_C)}{h}\\ \end{array} \right. \to \left\{ \begin{array}{c} \phi_A^e(\xi)=\frac{(x_C-x)(x_B-x)}{(x_C-x_A)(x_B-x_A)}\\ \phi_B^e(\xi)=\frac{(x_A-x)(x_B-x)}{(x_A-x_C)(x_B-x-C)}\\ \phi_C^e(\xi)=\frac{(x_A-x)(x_C-x)}{(x_A-x_B)(x_C-x_B)}\\ \end{array} \right.$$ Procedure: Compute $K_{kj}=\int_0^1 \bar{x}^{1/2}\phi_j'\bar{x}\phi_j'\bar{x}d\bar{x}$ after normalization $$\left[ \begin{array}{ccc} \int_0^h \sqrt{\bar{x}}\phi_1'(\bar{x})\phi_1'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_1'(\bar{x})\phi_2'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_1'(\bar{x})\phi_3'(\bar{x})d\bar{x} \\ \int_0^h \sqrt{\bar{x}}\phi_2'(\bar{x})\phi_1'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_2'(\bar{x})\phi_2'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_2'(\bar{x})\phi_3'(\bar{x})d\bar{x} \\ \int_0^h \sqrt{\bar{x}}\phi_3'(\bar{x})\phi_1'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_3'(\bar{x})\phi_2'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_3'(\bar{x})\phi_3'(\bar{x})d\bar{x} \\ \end{array} \right]$$ Compute $C_{kj}=2\int_0^1\phi_j(\bar{x})\phi_k(\bar{x})d\bar{x}$ $$ \left [ \begin{array}{ccc} \int_0^h \phi_1(\bar{x})\phi_1(\bar{x})d\bar{x} &\int_0^h \phi_1(\bar{x})\phi_2(\bar{x})d\bar{x} &\int_0^h \phi_1(\bar{x})\phi_3(\bar{x})d\bar{x} \\ \int_0^h \phi_2(\bar{x})\phi_1(\bar{x})d\bar{x} &\int_0^h \phi_2(\bar{x})\phi_2(\bar{x})d\bar{x} &\int_0^h \phi_2(\bar{x})\phi_3(\bar{x})d\bar{x} \\ \int_0^h \phi_3(\bar{x})\phi_3(\bar{x})d\bar{x} &\int_0^h \phi_3(\bar{x})\phi_2(\bar{x})d\bar{x} &\int_0^h \phi_3(\bar{x})\phi_3(\bar{x})d\bar{x} \\ \end{array} \right] $$ Solve for $\alpha_k$ in $(K_{kj}+C_{kj})\alpha_k = f_k \to \alpha_k=(K_{kj}+C_{kj})^{-1}f_k$ Approximated Solution: $\hat{u} \sim u_N=\sum_{k=1}^N \alpha_k(\bar{x})$

§ Discussion:

• Approximated solution using quadratic basic function converges to the analytical solution quicker than linear functions.
• We generally work with the weak form to minimize the residual
• The smaller the mesh, the more accurate $u(x,\phi)$ will be
• In the discussion of error analysis, quadratic algorithm yields a smaller error than linear due to a larger big O.

Torsional Analysis of Solid Elliptical Membrane

Suppose we have the schematic diagram below and the assoicated nodes.

We first note the symmetry of the membrane and simplify the model by considering only the first quadrant. Equation of Ellipse is $$(\frac{x}{a})^2+(\frac{y}{b})^2=1$$ We segment each of the quadrant into two individual meshes: Biquadratic Quadralateral and Bilinear Triangle
Rectangular Element: $Nodes [1 \;2 \;3 \;5 \;6 \;7 \;8 \;9 \;0]$
Triangular Element: $Nodes [3 \;4 \;12 \;11 \;10 \;7]$

§ Procedure

The following demonstrates the procedure in calculating the stiffness matrix and loading vector of the elements.
1 Choose $\Omega$ and $\Phi_j$, $j=1,2,...,N_e$ and specify the x-y coordinates $(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$ of nodal points of each element.
2 Specify a set of $N_i$ integration points $(\xi_l,\eta_l), \; l=1,2,...,N_l$ and quadrature weights for $\Omega$}
3 Calculate the values of $\Phi_j, \partial \Phi_j / \partial \xi, \partial \Phi_j / \partial \eta$ at the integration points.}
4 Calculate the values of $x=x(\xi,\eta), y=y(\xi,\eta)$ and their derivatives at the integration points.}
5 Calculate the values of the Jacobian and the functions $\partial \xi /\partial x, \partial \xi / \partial y, \partial \eta / \partial x, \partial \eta / \partial y$}
6 Compute $\partial \phi_j^e / \partial x$ and $\partial \phi_j^e / \partial y$}
7 Calculate the values of $k, b$ and $f$
We may skip this step if we assume $k=f=1$
8 Using the results of steps 3 to 7, calculate the values of the integrands at the integration points and multiply each by $w_i|\mathbf{J}(\xi_l,\eta_l)|$
9 Sum the numbers to obtain $k_{ij}^e$ and $f_i^e$

§ Numerical Integration

For the integration, we use numerical method to determine the weights of the Gaussian Quadrature.
$$\int_{-1}^1 f(x)dx=\sum_{i=1}^nw_if(t_i)$$ $$(n+1)L_{n+1}(t)-(2n+1)tL_n(t)+nL_{n-1}(t)=0, L_0(t)=1, L_1(t)=t$$ $$ L_{n+1}(t)=\frac{(2n+1)tL_n(t)-nL_{n-1}(t)}{n+1} $$ For all n values between 2 and 7, we have the following table

n	$\mathbf{L_{n+1}}$
2	$L_2(t)=\frac{3tL_1(t)-L_0(t)}{2}= \frac{3t^2-1}{2}$
3	$L_3(t)=\frac{5tL_1(t)-2L_1(t)}{3}= \frac{5t^3-3t}{2}$
4	$L_4(t)=\frac{7tL_1(t)-3L_2(t)}{4}= \frac{35t^4}{8}-\frac{15t^2}{4}+\frac{3}{8}$
5	$L_5(t)=\frac{63t^5}{8}-\frac{35t^3}{4}+\frac{15t}{8}$
6	$L_6(t)=\frac{231t^6}{16}-\frac{315t^4}{16}+\frac{105t^2}{16}-\frac{5}{16}$
7	$L_7(t)=\frac{429t^7}{16}-\frac{693t^5}{16}+\frac{315t^3}{16}-\frac{35t}{16}$

Graphical representation of the polynomials are below. Note that the x-intercepts are the roots of $L_n(t)$ over n.


Suppose $f(t)=t^k$ where $k=0,1,\ldots n$ $$\int_{-1}^1 f(t)dt = \int_{-1}^1 t^kdt = \frac{1-(-1)^{k+1}}{k+1} = \sum_{i=1}^n w_if(t_i)$$ Thus the integration weights can be determined by \begin{align}w_1t_1^k+\dots+w_nt_n^k=0 \qquad for \; k=1,3,\dots,2n-1\\ w_1t_1^k+\dots+w_nt_n^k=\frac{2}{k+1} \qquad for \; k=0,2,\dots,2n-2\end{align}

§ Torsional Properties

Physical quantities of interest, such as shear stresses and the relationship between the twisting moment, or torque, T, and the angle of twist $\theta$, per unit length of the shaft, can be determined as follows, once the potential function u is known. $$T =2G\theta \int_\Omega ud\Omega$$ The shear stresses on the elliptical cross-section are given by the following expressions: $$\sigma_{xz}=2G\theta \frac{\partial u}{\partial y}, \; \sigma_{yz}=-2G\theta \frac{\partial u}{\partial x}$$ The stress function may be written as $$\phi=B\lbrace(\frac{x}{a})^2+(\frac{y}{b})^2-1)\rbrace$$ But since $$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=-2G\theta$$ We get that $$B=-\frac{a^2b^2G\theta}{a^2+b^2}$$ $$ \sigma_{xz}=\frac{\partial \phi}{\partial y}=\frac{2By}{b^2}\quad \sigma_{yz}=-\frac{\partial \phi}{\partial x}=-\frac{2Bx}{a^2} $$ $$ T =2G\theta \int_\Omega ud\Omega=-\pi B ab $$ We can also verfiy that $$\phi=B\lbrace(\frac{x}{a})^2+(\frac{y}{b})^2-1)\rbrace=\phi=-\frac{a^2b^2G\theta}{a^2+b^2}\lbrace(\frac{x}{a})^2+(\frac{y}{b})^2-1)\rbrace$$

§ Procedure in Calculating Torsional Rigidity

$$u_h^e=\sum_{j=1}^{N_e}u_j^e \Phi_j^e$$ $$T =2G\theta \int_\Omega u_h^e d\Omega$$

1. With $u_h^e$ from part (i), multiply this with $\Phi_j$ for quadrilateral and triangular elements from page 198 and 204, respectively. You could ignore the nodes on boundary since it will be zero. For simplicity, consider working with one quadrant.
2. Use numerical integration to compute the integral similar to step 8 from part (i)
3. Repeat Steps 1 and 2 for the other element.
4. Sum everything from both elements.
5. Assume G and $\theta$ to be 1, multiply result from step (4) by 8 (2 from the constant term of the equation, and 4 since we have four quadrants.
6. Make sure to add node 1 three times since it's shared by all quadrants.

§ Discussion:

• Maximum shear stress occurs at the extreme values, namely at a and b (as it approaches the boundary). Furthermore, from Figure 9, we see that greater max stress occurs at the end of the minor axis of the ellipse
• Stress along the two axes of the centerlines is symmetrical due to the geometric symmetry of the membrane. This is expected since the heaviest concentration is the middle and distributed evenly across the membrane.

Stress Analysis on Rectangular Membrane

§ Derivation of Matrices

Evaluations of $K_{ij}$
$$K_{ij}^e=\int_{\Omega_e} k \left[\frac{\partial \phi_i}{\partial x} \frac{\partial \phi_y}{\partial x}+\frac{\partial \phi_i}{\partial y} \frac{\partial \phi_y}{\partial y} \right]$$ For rectangular elements $$\phi_1(x,y)=\frac{(x-a)(y-b)}{ab}$$ $$\phi_2(x,y)=\frac{-x(y-b)}{ab}$$ $$\phi_3(x,y)=\frac{xy}{ab}$$ $$\phi_4(x,y)=\frac{-(x-a)y}{ab}$$ Evaluations of $f_i$
$$f_i^e=\int_{\Omega_e} f\phi_idxdy=\frac{1}{2A_e}\int_{\Omega_e}f(x,y)(\alpha_i +\beta_ix+\gamma_iy)dxdy$$ From the Law of Conservation Energy, the strain and potential energy must be equal and opposite from each other
$S=-W=\frac{1}{2}u^{h^T} ku^h$

§ Sample Output

For a fine mesh, students should see the greater yield of higher stress distribution from the centroid of the membrane.
Hint: Play with the mesh size on the element slider of the GUI. What would happen if there was a slit (cut) on the membrane?

§ Discussion:

• The 4 by 4 matrix elementary stiffness matrix is the same for any given n.
• Stress deflection along the two axes of the centerlines is symmetrical due to the geometric symmetry of the membrane
• Introduction of a slit on any plane causes an uneven symmetry about both axes
• Strain and potential energy derivation is consistent with the law of conservation of energy. Membrane is able to “store” more energy as we refine better meshes
• Suppose a slit cut is on the left of the membrane, because the distribution of energy originates from the center of the membrane, one can presume that the stress on the right is always greater than the left.

Thermal Analysis on Triangular Membrane

Steady-state Heat Transfer Using Tri-Quadratic Hexahedral Finite Element

The 27-node hexahedron is the analog of the 8-node “serendipity” quadrilateral
For example, the general formulas for the midside nodes are $$N_j=\frac{1}{4}(1-\xi^2)(1+\eta_j\eta)(1+\varsigma_j\varsigma)$$ $$N_j=\frac{1}{4}(1+\xi_j\xi)(1-\eta^2)(1+\varsigma_j\varsigma)$$ $$N_j=\frac{1}{4}(1+\xi_j\xi)(1+\eta_j\eta)(1-\varsigma^2)$$ Hexahedral element with tri-quadratic approximation functions
$N_1=\frac{1}{8}(1-\xi)(1-\eta)(1-\varsigma)$
$N_2=\frac{1}{8}(1+\xi)(1-\eta)(1-\varsigma)$
$N_3=\frac{1}{8}(1+\xi)(1+\eta)(1-\varsigma)$
$N_4=\frac{1}{8}(1-\xi)(1+\eta)(1-\varsigma$
$N_5=\frac{1}{8}(1-\xi)(1-\eta)(1+\varsigma)$
$N_6=\frac{1}{8}(1+\xi)(1-\eta)(1+\varsigma)$
$N_7=\frac{1}{8}(1+\xi)(1+\eta)(1+\varsigma)$
$N_8=\frac{1}{8}(1-\xi)(1+\eta)(1+\varsigma)$
$N_9=\frac{1}{4}(1-\xi^2)(1-\eta)(1-\varsigma)$
$N_{10}=\frac{1}{4}(1+\xi)(1-\eta^2)(1-\varsigma)$
$N_{11}=\frac{1}{4}(1-\xi^2)(1+\eta)(1-\varsigma)$
$N_{12}=\frac{1}{4}(1-\xi)(1-\eta^2)(1-\varsigma)$
$N_{13}=\frac{1}{4}(1+\xi)(1-\eta)(1+\varsigma^2)$
$N_{14}=\frac{1}{4}(1+\xi)(1-\eta^2)(1+\varsigma^2)$
$N_{15}=\frac{1}{4}(1-\xi^2)(1+\eta)(1+\varsigma)$
$N_{16}=\frac{1}{4}(1-\xi)(1-\eta^2)(1+\varsigma)$
$N_{17}=\frac{1}{4}(1-\xi)(1-\eta)(1-\varsigma^2)$
$N_{18}=\frac{1}{4}(1+\xi)(1-\eta^2)(1-\varsigma^2)$
$N_{19}=\frac{1}{4}(1+\xi)(1+\eta)(1-\varsigma^2)$
$N_{20}=\frac{1}{4}(1-\xi)(1+\eta^2)(1-\varsigma^2)$
$N_{21}=\frac{1}{2}(1-\xi^2)(1-\eta^2)(1-\varsigma)$
$N_{22}=\frac{1}{2}(1-\xi^2)(1-\eta^2)(1+\varsigma)$
$N_{23}=\frac{1}{2}(1-\xi^2)(1-\eta)(1-\varsigma^2)$
$N_{24}=\frac{1}{2}(1+\xi^2)(1-\eta^2)(1-\varsigma^2)$
$N_{25}=\frac{1}{2}(1-\xi^2)(1+\eta)(1-\varsigma^2)$
$N_{26}=\frac{1}{2}(1-\xi)(1-\eta^2)(1-\varsigma^2)$
$N_{27}=(1-\xi^2)(1-\eta^2)(1-\varsigma^2)$

Derivatives of Shape Functions \begin{eqnarray} \frac{\partial N_i^{(e)}}{\partial x}=\frac{\partial N_i^{(e)}}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial N_i^{(e)}}{\partial \eta}\frac{\partial \eta}{\partial x}+\frac{\partial N_i^{(e)}}{\partial \varsigma}\frac{\partial \varsigma}{\partial x}\\ \frac{\partial N_i^{(e)}}{\partial y}=\frac{\partial N_i^{(e)}}{\partial \xi}\frac{\partial \xi}{\partial y}+\frac{\partial N_i^{(e)}}{\partial \eta}\frac{\partial \eta}{\partial y}+\frac{\partial N_i^{(e)}}{\partial \varsigma}\frac{\partial \varsigma}{\partial y}\\ \frac{\partial N_i^{(e)}}{\partial z}=\frac{\partial N_i^{(e)}}{\partial \xi}\frac{\partial \xi}{\partial z}+\frac{\partial N_i^{(e)}}{\partial \eta}\frac{\partial \eta}{\partial z}+\frac{\partial N_i^{(e)}}{\partial \varsigma}\frac{\partial \varsigma}{\partial z} \end{eqnarray} The infinitesimals $d\xi$, $d\eta$ and $d\varsigma$ transform into $dx$, $dy$, and $dz$ can be written in matrix form as $$ \left[ \begin{array}{c} dx \\ dy \\ dz \end{array} \right]= \left[ \begin{array}{ccc} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} & \frac{\partial x}{\partial \varsigma}\\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} & \frac{\partial y}{\partial \varsigma}\\ \frac{\partial z}{\partial \xi} & \frac{\partial z}{\partial \eta} & \frac{\partial z}{\partial \varsigma} \end{array} \right] \left[ \begin{array}{c} d\xi \\ d\eta \\ d\varsigma \end{array} \right] $$ Jacobian Matrix $$ \mathbf{J}=\frac{\partial (x,y,z)}{\partial (\xi,\eta,\varsigma)}= \left[ \begin{array}{ccc} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi}\\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta}\\ \frac{\partial x}{\partial \varsigma} & \frac{\partial y}{\partial \varsigma} & \frac{\partial z}{\partial \varsigma} \end{array} \right] $$ The isoparametric definition of hexahedron element geometry is $$ x=\sum_{i=1}^{20}x_iN_i^{(e)} \qquad y=\sum_{i=1}^{20}y_iN_i^{(e)} \qquad z=\sum_{i=1}^{20}z_iN_i^{(e)} $$ Whenever $|J|\neq 0$, we can write $$ \left[ \begin{array}{ccc} \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} & \frac{\partial \xi}{\partial z}\\ \frac{\partial \eta}{\partial x} & \frac{\partial \eta}{\partial y} & \frac{\partial \eta}{\partial z}\\ \frac{\partial \varsigma}{\partial x} & \frac{\partial \varsigma}{\partial y} & \frac{\partial \varsigma}{\partial z}\\ \end{array} \right] = \mathbf{J^{-1}} $$ The derivatives of $\Phi_j^e$ are obtained by the chain rule: \begin{eqnarray} \frac{\partial \Phi_j^e}{\partial x}=\frac{\Phi_j}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\Phi_j}{\partial \eta} \frac{\partial \eta}{\partial x}+\frac{\Phi_j}{\partial \varsigma} \frac{\partial \varsigma}{\partial x}\\ \frac{\partial \Phi_j^e}{\partial y}=\frac{\Phi_j}{\partial \xi} \frac{\partial \xi}{\partial y}+\frac{\Phi_j}{\partial \eta} \frac{\partial \eta}{\partial y}+\frac{\Phi_j}{\partial \varsigma} \frac{\partial \varsigma}{\partial y}\\ \frac{\partial \Phi_j^e}{\partial z}=\frac{\Phi_j}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\Phi_j}{\partial \eta} \frac{\partial \eta}{\partial z}+\frac{\Phi_j}{\partial \varsigma} \frac{\partial \varsigma}{\partial z}\\ \end{eqnarray} Numerical Integration Over Hexahedral, the total number of Gauss points is $p^3$ $$ \int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}F(\xi,\eta,\varsigma)d\xi d\eta d\varsigma \approx \sum_{i=1}^{p_1}\sum_{j=1}^{p_2}\sum_{k=1}^{p_3} w_iw_jw_k F(\xi_i,\eta_j,\varsigma_k) $$

§ Bending of a Uniform, Homogeneous Elastic Beam (Euler-Bernoulli Theory)

Strong Form: $$\int_0^L\lbrace (EIw'')''-q\rbrace vdx-\lbrace -(EIw'')(L)-M_L\rbrace v'(L)+\lbrace (EIw'')(L)-v_L\rbrace v(L)=0 \quad \forall \;v \quad s.t \; v(0)=0, \; v'(0)=0$$ Weak Form: $$\int_0^L EIw''v''dx-\int_0^L qvdx+M_Lv'(L)-v_Lv(L)=0 \quad \forall \; v \quad s.t \; v(0)=0, \; v'(0)=0$$ Principal of Virtual Work $$\int_0^L EIw''v''dx=\int_0^L qvdx-M_Lv'(L)+v_Lv(L)$$ But since we don't have any applied moment of load at the supported end,our governing equation is actually $$\int_0^L EIw''v''dx=\int_0^L qvdx$$ The simplest Bernoulli-Euler plane beam element with two end nodes has four degrees of freedom $\mathbf{u^e}=[v_1 \theta_1 v_2 \theta_2]^T$.
Shape functions for this problem are conveniently expressed in terms of the dimensionless coordinate $$\xi=\frac{2x}{h}-1 \qquad \frac{dx}{d\xi}=\frac{1}{2}h \qquad \frac{d\xi}{dx}=\frac{2}{h}$$ \begin{eqnarray} N_1(\xi)=\frac{1}{4}(1-\xi)^2(2+\xi)\\ N_2(\xi)=\frac{h}{8}(1-\xi)^2(1+\xi)\\ N_3(\xi)=\frac{1}{4}(1+\xi)^2(2+\xi)\\ N_4(\xi)=\frac{h}{8}(1+\xi)^2(\xi-1)\\ \end{eqnarray} Local Stiffness Matrix $$K_{ij}=\int_eEI\frac{d^2N_i(x)}{dx^2} \frac{d^2N_j(x)}{dx^2}dx=\int_{-1}^1 EI\frac{d^2N_i(x)}{dx^2} \frac{d^2N_j(x)}{dx^2}\frac{1}{2}h d\xi$$ For a generic case with both free nodes, the matrix becomes \[ K^e = \frac{EI}{h^3}\left[ \begin{array}{cccc} 12 & 6h & -12 & 6h \\ 6h & 4h^2 & -6h & 2h^2 \\ -12 & -6h & 12 & -6h \\ 6h & 2h^2 & -6h & 4h^2 \end{array} \right]\] Local Force Vector $$f_i^{(e)}=\int_e q(x)N_i(x)dx=\int_{-1}^1 q(x)N_i(x)\frac{1}{2}h d\xi$$ For uniform load $q_0$ $$f^{(e)}=\frac{q_0h_e}{12} \left\{ \begin{array}{c} 6\\ h_e \\ 6 \\ -h_e \end{array} \right\} $$ System array is defined with dimensions $n\;x\;m$ where $n$ is 4 (the number of degree of freedom per element) and $m$ is 3 (the number of elements).

§ Distribution of Bending Moment Along the Beam

Since moment is related to the distributed load by its second derivative, the approximate solution is $$M(x)=-EI\frac{d^2w}{dx^2}=-EI\sum_{j=1}^4 u_j^e \frac{d^2 \Phi_j^e}{dx^2}$$

§ Reactions at the Beam Ends

$$V(x)=\frac{dM}{dx}=-\frac{d}{dx}\left(EI\frac{d^2w}{dx^2}\right)=-EI\sum_{i=1}^4 u_j^e \frac{d^3 \Phi_j^e}{dx^3}$$

1. [C#] Given a string, find the most frequent character
2. [C#] Reverse a given string
3. [Java] Anagrams
4. [Java] Array manipulation
5. [Python] Majority element
6. Java] Detect whether a linked list contains-a-cycle
7. [Python] Mini max sum
8. [Java] Grading students
9. [Java] Beautiful days at the movies
10. [Java] Palindrome
11. [C#] Recursive digit sum
12. [Java] Sparse arrays
13. [Java] Lexicographical substring comparison
14. [Python] Move the spaces
15. [Java] Two sum
16. [C#] Caesar Cipher
17. [C#] Max min sum
18. [C#] Power sum
19. [Java] Cut the sticks
20. [Java] Sock merchants
21. [Python] Lily's homework
22. [Java] Diagonal difference
23. [C#] Two sum for sorted array
24. [Java] Viral advertising
25. [C#] Plus one
26. [Java] Substring comparisons
27. [C#] Same tree

Windows Presentation Foundation (or WPF) is a graphical subsystem by Microsoft for rendering user interfaces in Windows-based applications. In my coding career I worked with this technology. This page provides a mini project that has a few comprehensive examples of what WPF is capable of. The last topic of this page talks about Model-View-ViewModel (MMVV), which is a software architectural patttern that decouples the user interface with the business logic. This page is not meant to be a tutorial but simply a refresher to those that have experiences with WPF and MVVM.

Click on the icon to download the Visual Studio solution of this demo project

This project demostrates a few WPF features through a series of examples. The first example is on animation. Here I have two images, one is a normal picture of me and the another one is not so skinny picture of me. These two images are overlapped with each other using the Canvas tag. In WPF, every control can only have one child, thus containers like StackPanel, WrapPanel, and DockPanel are used to lay out the interface, howevever they are "stacked" relative to each other except for the canvas and grid containers. Only Grid and Canvas allow for overlapping. Here we use a trigger with a loaded event to begin an animation that changes the opacity of a target with the name "pic" (which in this case is the not skinny picture) from fully visible to invisible within two seconds. The "RepeatBehavior" is set to forever which means the animation loops infinitely. Furthermore, we have two more routed events for mouse in and out which pauses and resumes the storyboard, respectively.

                            
<Canvas>
    <Canvas.Triggers>
        <EventTrigger RoutedEvent="Canvas.Loaded">
            <BeginStoryboard x:Name="OpacityStoryboard">
                <Storyboard TargetName="pic2" TargetProperty="Opacity">
                    <DoubleAnimation From="1" To="0" Duration="0:0:2" AutoReverse="True" RepeatBehavior="Forever" />
                </Storyboard>
            </BeginStoryboard>
        </EventTrigger>

        <EventTrigger RoutedEvent="UIElement.MouseEnter">
            <PauseStoryboard BeginStoryboardName="OpacityStoryboard" />
        </EventTrigger>
        <EventTrigger RoutedEvent="UIElement.MouseLeave">
            <ResumeStoryboard BeginStoryboardName="OpacityStoryboard" />
        </EventTrigger>
    </Canvas.Triggers>

    <Image x:Name="pic1" Source="Images/Skinny.jpg" />
    <Image x:Name="pic2" Source="Images/NotSkinny.jpg" />
</Canvas>
                            
                            

In this demo, we observe some pretty amazing WPF styling. The first style targets the whole Window. It has a trigger that listens to the binding of the property "IsChecked" of the element by the name of "redColorCheckBox." If this value is true, then it sets the background to red. But what background? The background of the specific type, which in this case is the window. The second block of code is inside the Window.Resources which means any styling or defines in this block is applied to the whole window. The first style in here applies to every button, hence in the video clip on the left you see that the three buttons all look the same (height of 30px, height of 80px, font size of 12px, etc). Furthermore, once again we have a trigger and this one is "IsMouseOver", meaning whenever the mouse is over the control (button), the foreground turns red. Styles in WPF do not have to be defined under Resources, it can be defined straight within the scope of the individual control. If you look at the second textbox, you will see that I am manipulating the "Background" and "IsEnabled" properties. The background is binded to whatever value I type in the textbox. This is done through RelativeSource to myself and taking the text property. Similarly, I have a data trigger that compares to my typed text, if it is equal to the string "disabled" then it will set the "IsEnabled" property to false. As you can see, the possibilities of WPF styling is endless. The complete code snippet is below.

  
<Window.Style>
	<Style TargetType="{x:Type Window}">
		<Style.Triggers>
			<DataTrigger Binding="{Binding ElementName=redColorCheckBox, Path=IsChecked}" Value="True">
				<Setter Property="Background" Value="Red" />
			</DataTrigger>
		</Style.Triggers>
	</Style>
</Window.Style>

<Window.Resources>
	<Style TargetType = "{x:Type Button}">
		<Setter Property = "Height" Value = "30" />
		<Setter Property = "Width" Value = "80" />
		<Setter Property = "Foreground" Value = "Blue" />
		<Setter Property = "FontSize" Value = "12" />
		<Setter Property = "Margin" Value = "10" />
		<Style.Triggers>
			<Trigger Property="IsMouseOver" Value="True">
				<Setter Property="Foreground" Value="red" />
			</Trigger>
		</Style.Triggers>
	</Style>

	<Style x:Key="practiceStyle" TargetType="{x:Type TextBox}" >
		<Style.Triggers>
			<MultiTrigger>
				<MultiTrigger.Conditions>
					<Condition Property="IsMouseOver" Value="True" />
					<Condition Property="IsKeyboardFocused"  Value="True" />
				</MultiTrigger.Conditions>
				<MultiTrigger.Setters>
					<Setter Property="Foreground" Value="red" />
					<Setter Property="FontWeight" Value="Bold" />
					<Setter Property="FontSize" Value="18" />
				</MultiTrigger.Setters>
			</MultiTrigger>
		</Style.Triggers>
	</Style>

	<SolidColorBrush x:Key="staticBrush" Color="Yellow" />
</Window.Resources>

<StackPanel HorizontalAlignment="Center">
	<Button Content = "Button1"/>
	<Button Content = "Button2"/>
	<Button Content = "Button3"/>

	<TextBox Width="150" Background="{StaticResource staticBrush}" Style="{StaticResource practiceStyle}" />

	<TextBox Width="100" Margin="10">
		<TextBox.Style>
			<Style TargetType="TextBox">
				<Style.Triggers>
					<DataTrigger Binding="{Binding RelativeSource={RelativeSource Self}, Path=Text}" Value="disabled">
						<Setter Property="IsEnabled" Value="False" />
					</DataTrigger>
				</Style.Triggers>

				<Setter Property="Background" Value="{Binding RelativeSource={RelativeSource Self}, Path=Text}" />
			</Style>
		</TextBox.Style>
	</TextBox>

	<CheckBox x:Name = "redColorCheckBox" Content = "Set red as foreground color" Margin = "10"/>
</StackPanel>
  

Data binding is what makes WPF so powerful. In this demo, we will see how to properly create a data template. The XAML is pretty straight forward, we have a listbox and a button. This button subscribes to a click event which we will discuss later. The listbox contains an ItemsSource that populates the collection as well as an Item Template. The StaticResource of this property points to a data template that "customizes" the appearance on how to display the data, which in this sense I assign it to have a grid inside a border and in this grid we get and set the name and age of a person. Note that the binding of "Name" and "Age" are actually properties of an instance of the datum (from the item source). As for the population of the items and the model itself, I handled them in the code behind. I have a class called Person that has name and age as attributes and publicly expose them so the data template can retrieve them. In the constructor I have a list of people and assign it to the item source of the listbox. Note that we can also use the .NET ObservableCollection which implements the INotifyCollectionChanged that works well with MVVM. When users change the name and age of each item, the properties automatically get updated because of the binding, and when they click on the button, there's a message box that shows the person's information.

                            
<Window.Resources>
	<DataTemplate x:Key ="template">
		<Border x:Name="bord3r" BorderBrush="Red" BorderThickness="1">
			<Grid>
				<Grid.RowDefinitions>
					<RowDefinition Height = "Auto" />
					<RowDefinition Height = "Auto" />
				</Grid.RowDefinitions>
				<Grid.ColumnDefinitions>
					<ColumnDefinition Width = "Auto" />
					<ColumnDefinition Width = "200" />
				</Grid.ColumnDefinitions>

				<Label Margin = "10" Content="Name"/>
				<TextBox Grid.Column = "1" Margin = "10" Text = "{Binding Name}" />
				<Label Margin = "10" Grid.Row = "1" Content="Age"/>
				<TextBox Grid.Column = "1" Grid.Row = "1" Margin = "10" Text = "{Binding Age}" />
			</Grid>
		</Border>

		<DataTemplate.Triggers>
			<DataTrigger Binding="{Binding Path=Name}" Value="Yvonne">
				<Setter TargetName="bord3r" Property="BorderBrush" Value="blue" />
			</DataTrigger>
		</DataTemplate.Triggers>
	</DataTemplate>
</Window.Resources>

<Grid>
	<Grid.RowDefinitions>
		<RowDefinition Height = "Auto" />
		<RowDefinition Height = "*" />
	</Grid.RowDefinitions>

	<ListBox x:Name="listbox" ItemsSource = "{Binding Source}" ItemTemplate="{StaticResource template}" />

	<Button Grid.Row = "1"  Content = "_Show..." Click = "Button_Click" Width = "80" HorizontalAlignment = "Left" Margin = "10"/>

</Grid>
    
    
public partial class DataTemplateExample : Window
{
    public DataTemplateExample()
    {
        InitializeComponent();

        List people = new List();
        people.Add(new Person { Name = "Ethan", Age = 27 });
        people.Add(new Person { Name = "Yvonne", Age = 62 });
        people.Add(new Person { Name = "Thomas", Age = 12 });

        listbox.ItemsSource = people;
    }

    private void Button_Click(object sender, RoutedEventArgs e)
    {
        Person selectedPerson = (Person)listbox.SelectedValue;

        if (selectedPerson != null)
        {
            string message = string.Format("{0} is {1} years old", selectedPerson.Name, selectedPerson.Age);
            MessageBox.Show(message);
        }
    }
}

public class Person
{
    private string _Name;
    private double _Age;

    public string Name
    {
        get { return _Name; }
        set { _Name = value; }
    }

    public double Age
    {
        get { return _Age; }
        set { _Age = value; }
    }
}
                            
                            

Now this is a real treat. In 2016 there was a big PowerBall lottery of a half a billion jackpot. I was inspired by the event to write a simple app just for fun. This side project utilizes many WPF and .NET features ranging from storyboard animation to multi-threading. Click on the icon to download the Visual Studio solution of this Powerball project

This is an example of a very basic and classic MVVM pattern. The main idea of MVVM is that the Model should know nothing about the View and vice-versa. The Model is defined as any object that holds information. The View is the front end presentational layer. The "link" between the two is the ViewModel. The ViewModel should only know about the Model and not the View, and the View should only know about the ViewModel and not the Model.

This example is very simple, we have a listbox, data grid, and combo box that shares the same items source. There is a button that upon invoke will take a textbox string and adds to the collection. The core of MVVM lies in the implementation of the INotifyPropertyChanged interface. This interface allows any messages to be updated back to the View. Any property in the ViewModel that is bound to the View should implement this.

                            
public class ViewModelBase: INotifyPropertyChanged
{
    public event PropertyChangedEventHandler PropertyChanged;

    protected virtual void OnPropertyChanged(string propertyName)
    {
        PropertyChangedEventHandler handler = PropertyChanged;
        if (handler != null)
        {
            handler(this, new PropertyChangedEventArgs(propertyName));
        }
    }
}
                            
                            

The second main component is the implementation of the ICommand interface. This is to bind commands in the View such as button or any control event.

                            
public class DelegateCommand : ICommand
{
    private readonly Action _execute;
    private readonly Func _canExecute;

    public DelegateCommand(Action executeMethod)
        : this(executeMethod, null)
    {
    }

    public DelegateCommand(Action executeMethod, Func canExecuteMethod)
    {
        if (executeMethod == null)
            throw new ArgumentNullException("executeMethod");

        _execute = executeMethod;
        _canExecute = canExecuteMethod;
    }

    public void Execute(object parameter)
    {
        _execute();
    }

    public bool CanExecute(object parameter)
    {
        return _canExecute == null ? true : _canExecute();
    }

    public event EventHandler CanExecuteChanged
    {
        add { CommandManager.RequerySuggested += value; }
        remove { CommandManager.RequerySuggested -= value; }
    }
}
                            
                            

The model can be anything, here I choose to create a blueprint of a person that has a first name, last name, and age.

                            
public class Person : INotifyPropertyChanged
{
    private string _FirstName;
    private string _LastName;
    private int _Age;

    public string FirstName
    {
        get { return _FirstName; }
        set
        {
            if (_FirstName != value)
            {
                _FirstName = value;
                OnPropertyChanged("FirstName");
            }
        }
    }

    public string LastName
    {
        get { return _LastName;}
        set
        {
            if (_LastName != value)
            {
                _LastName = value;
                OnPropertyChanged("LastName");
            }
        }
    }

    public int Age
    {
        get { return _Age;}
        set
        {
            if (_Age != value)
            {
                _Age = value;
                OnPropertyChanged("Age");
            }
        }
    }

    public event PropertyChangedEventHandler PropertyChanged;

    protected virtual void OnPropertyChanged(string propertyName)
    {
        PropertyChangedEventHandler handler = PropertyChanged;
        if (handler != null)
        {
            handler(this, new PropertyChangedEventArgs(propertyName));
        }
    }
}
                            
                            

I daresay the ViewModel is the most complex part of the application. Afterall, it handles all the business logic and serves as the mediator between the Model and View. Here, we populate the collection in the constructor and listen to the change in selected person. Because the list box, data grid, and combo box all shared the same ObservableCollection and because we notify the same messages back to the view, all three get updated at the same time.

                            
public class MainWindowViewModel : ViewModelBase
{
    public DelegateCommand AddUserCommand { get; set; }
    public ObservableCollection People { get; set; }
    private Person _SelectedPerson;

    private string _SelectedItemString;
    public string TextProperty { get; set; }

    public MainWindowViewModel()
    {
        AddUserCommand = new DelegateCommand(OnAddUserCommand);

        People = new ObservableCollection
        {
            new Person { FirstName="Ethan", LastName="Uong", Age=32 },
            new Person { FirstName="Yvonne", LastName="Liu", Age=26 },
            new Person { FirstName="Happy", LastName="Doggy", Age=3 },
        };
    }

    public Person SelectedPerson
    {
        get { return _SelectedPerson; }
        set
        {
            if (_SelectedPerson != value && value != null)
            {
                _SelectedPerson = value;
                SelectedItemString = value.FirstName;
                OnPropertyChanged("SelectedPerson");
            }
        }
    }

    public string SelectedItemString
    {
        get { return _SelectedItemString; }
        set
        {
            if (_SelectedItemString != value)
            {
                _SelectedItemString = value;
                OnPropertyChanged("SelectedItemString");
            }
        }
    }

    private void OnAddUserCommand()
    {
        if (!string.IsNullOrEmpty(TextProperty))
        {
            People.Add(new Person {
                FirstName = TextProperty.ToString(),
                LastName = TextProperty.ToString(),
                Age = DateTime.Now.Second
            });
        }
    }
}
                            
                            

And finally we have the View. The last major concept is how exactlyl do we bind everything together? We have established that the View must know nothing about the Model, but how precisely does the View understand the ViewModel? This is done by setting the DataContext of the View to an instance of the ViewModel.

                            
<Grid Margin="20">
    <Grid.RowDefinitions>
        <RowDefinition Height="Auto"/>
        <RowDefinition Height="Auto"/>
    </Grid.RowDefinitions>

    <StackPanel Grid.Row="0">
        <StackPanel Orientation="Horizontal">
            <ListBox ItemsSource="{Binding People}" SelectedItem="{Binding SelectedPerson}"
                     DisplayMemberPath="FirstName" HorizontalAlignment="Left"/>
            <DataGrid ItemsSource="{Binding People}" SelectedItem="{Binding SelectedPerson}" CanUserAddRows="False"
                      HorizontalAlignment="Left" Margin="5,0,0,0"/>
            <ComboBox ItemsSource="{Binding People}" SelectedItem="{Binding SelectedPerson}"
                      DisplayMemberPath="FirstName" Margin="5,0,0,5" VerticalAlignment="Top"/>
        </StackPanel>

        <TextBlock FontWeight="Bold" Margin="5" Text="The selected person is ">
        <Run Text="{Binding SelectedItemString}"/></TextBlock>
        <Label Content="Type in a name and hit button to add to collection" />
    </StackPanel>

    <StackPanel Grid.Row="1" Width="150" HorizontalAlignment="Left">
        <TextBox Text="{Binding TextProperty}" Margin="5"/>
        <Button Content="Add person" Command="{Binding AddUserCommand}" Margin="5" />
    </StackPanel>
</Grid>
                            
                            

Click on the icon to download the Visual Studio solution of this basic MVVM example

Abstract

The study of Chaos Theory to analyze unprecedented events has been advancing since mid 20th century. Scientists and researchers find many applications from this field such as weather rediction and explanation of the rise and fall of stocks through the examination of the Lorenz attractor. This page seeks to have a better understanding of the chaotic behavior of the deterministic system by solving the system of ordinary differential equations using two numerical methods (Runge Kutta to the 4th order and the adaptive Runge Kutta. It is shown that the adaptive method yields significant better results and the analysis of the Lorenz attractor plays a vital role in understanding chaos theory. For my numerical method grad course, I designed a software package to understand the behavior of the Lorenz attractor on different parameters.

Introduction

Chaos theory is a field of study in mathematics that concerns the behavior of dynamical systems with sensitive initial conditions. It is based on the notion that a small difference in initial conditions such as loss of significance and rounding errors can cause diverging outcomes to the system (deterministic or stochastic), thus rendering it predictable in the long term. Chaos, as summarized by Edward Lorenz, is "when the present determines the future, but the approximate present does not approximately determine the future [1]." Chaos theory deals with nonlinear systems that are impossible to predict or control such as turbulence, weather, brain states, and stock market. One main characteristic of this theory is its unpredictability. Because we can never know if the initial conditions are sufficient enough, the ultimate fate of the system cannot be determined. A slight error in measuring the state of the system can dramatically alter the result, and as such long-range weather prediction is not possible. Chaos exists in many natural phenomena such as weather and one approach to mathematically model its behavior is through the Lorenz attractor. The overall idea of chaos theory can be summarized by Poincar´e in 1903: "A very small cause which escapes our notice determines a considerable effect that we cannot fail to see...even if the case that the natural laws had no longer secret for us...we could only know the initial situation approximately...It may happen that small differences in initial conditions produce very great ones in the final phenomena."

Motivation

The motivation in this study is that chaotic systems have been widely adopted in many areas of study such as chemical reaction, weather forecast, electrical circuit, stockmarket, human biology , message encryption, etc. The subsequent paragraphs in this motivation section is devoted to explaining an example of the application of chaotic systems in message encryption. One concept that one needs to grasp in order to understand how chaotic systems work in message encryption is that in the chaotic world it is impossible to build two chaotic systems to produce the same output. However, ”if two identical stable systems are driven by the same chaotic signal,” the outputs are two chaotic signals that are identical to each other. With this concept in mind, the message encryption process can be presented as the following [3].

Referring to the figure below, at the sending location chaotic signal x is generated and fed into a stable system F to produce the output as another chaotic signal y. Message encryption occurs when the message is mixed with the chaotic signal y and the result is the message-carrying signal y* that has a chaotic behaviour. Next step is to transmit two signals to the receiving location. One is the original chaotic signal x and the other is the message-carrying chaotic signal y*. At the sending location only chaotic signal x is fed into another stable system F which is identical to the stable system F at the sending location. Therefore, the output will be the chaotic signal y, which is the same output produced by the corresponding stable system F at the sending location. Finally, the original message can be recovered by subtracting chaotic signal y from the message-carrying chaotic signal y* at the receiving location. [3]

Interested Concentration

A Lorenz system is a system of ordinary differential equations with chaotic solutions for certain parameter values and initial conditions. When plotted, the system resembles a butterfly, hence chaos theory is often time called the Butterfly Effect. A simple analogy is the following: a flap of a butterfly in Mexico can cause a tornado in Texas. The notion is that a slight change in the current state could potentially cause a drastic change in the long-term effect, thus the validity of the starting points is highly sensitive.

In Lorenz system, $x$, $y$, and $z$ are the system state, $t$ is time, $\sigma$, $\rho$, and $\beta$ are the system parameters. The constants $\rho$, $\beta$, and $\sigma$ determine the behavior of the system. Lorenz defines a deterministic sequence as one in which only one thing can happen next from an existing event. Furthermore, deterministic chaos is something that looks random, but is still deterministic. For the purpose of this page, we will use the canonical parameters $\rho=10$, $\beta=\frac{8}{3}$, and $\sigma=28$.

Mathematics and Procedure

Any $n^{th}$ order linear ordinary differential equation can be expressed in the form of a system of first order linear ODEs. Numerical approximation is applied to ascertain values that may otherwise be impossible to solve through traditional methods. One popular family of numerical methods is the Runge Kutta which carries a better computational accuracy than the other methods [4]. Runge Kutta method consists of a number of numerical methods that can provide different levels of accuracy. Runge Kutta Fourth Order is one of the Runge Kutta methods and was introduced in the class. Therefore, the goal of this study is to introduce a new Runge Kutta method called Adaptive Runge ­Kutta method that is hypothetically capable of providing a better accuracy and efficiency than Runge ­Kutta Fourth Order method. The procedure in this study is to solve the Lorenz attractor by using both Adaptive Runge ­Kutta and RK4 methods and compare the results. The motivation for comparing the two methods is that the efficiency gained in adaptive RK can be hundred times or more compared to RK4 [5].

Nomenclature

Runge Kutta Fourth Order

As presented in the class, Runge ­Kutta Fourth Order or RK4 is one amongst a family of numerical methods developed by two German Mathematicians C. Runge and M.W. Kutta that are cable of approximating the solutions to ordinary differential equations by discretizing both spatial and temporal domains. Runge ­Kutta methods can also be extended to solve a system of initial value problems (IVPs) ­ which is the case for the Lorenz Attractor. RK4 requires four ­function evaluations as presented in equation 4 to 8 and it has an error in the order of $O(h^4)$ [6].

Adaptive Runge Kutta

Adaptive Runge ­Kutta method is the extension of RK4 method in that RK4 uses a uniform time step size throughout the entire domain (Figure 2.a) while Adaptive Runge­Kutta method uses non-uniform time step size by employing the step size control method in order to obtain an estimate within a prescribed tolerance $\epsilon$ and with a minimum computational effort. The purpose of adapting the step size is to use smaller steps with steep slope regions of the solution curve and to use larger steps with shallow slope regions of the solution curve in order to obtain an estimate within $\epsilon$

In order to adapt the step size, expected error estimate is needed to compare to the prescribed tolerance, $\epsilon$. If the expected error is smaller or equal to the tolerance, the current step size shall be used for the next step otherwise a new step size is to be computed based on the estimate error. There are two common methods in estimating the expected error. The first method, called step doubling, is to take each step twice using step sizes $h$ as full step and $h/2$ as two small steps. The expected error is calculated as the difference between the two results. The second method is to use higher-­lower order method in which the idea is to calculate the error as the difference between results from two order methods, say 4th and 5th, and adapt a unique step size for the next time step throughout.The adaptive Runge ­Kutta method conducted in this study is the Runge ­Kutta ­Fehlberg method, also known as Embedded Runge ­Kutta Formulas (RK45) [5].

Runge ­Kutta ­Fehlberg Method (RK45) was developed by a German mathematician Erwin Fehlberg in 1969 with error in order of $O(h^5)$. RK45 requires six function evaluations. A combination of these functions provides the calculation of the next step based on RK4, $x_{i+1}^{RK4}$, (eq.15) and a different combination of these functions provides the calculation of the next step based on RK5, $x_{i+1}^{RK5}$

If $R <=\epsilon$ keep $x$ as the current step solution and move to the next step with step size $\delta h$.

If $R>\epsilon$ recalculate the current step with step size $\delta h$.

Step-doubling as a means for adaptive stepsize control in fourth-order Runge-Kutta. Points where the derivative is evaluated are shown as filled circles. The open circle represents the same derivatives as the filled circle immediately above it, so the total number of evaluations is 11 per two steps. Comparing the accuracy of the big step with the two small steps gives a criterion for adjusting the stepsize on the next step, or for rejecting the current step as inaccurate [5].

Result and discussion

The solutions to Lorenz Attractor found using both methods are presented in Figure 4. These results are projected onto three principal planes xy, yz and xz for the purpose of comparison. It can be observed that in general the plots of solution after a given time, t, from RK45 are smoother than those from RK4. This indicates that RK45 is capable in providing a more accurate result than RK4.

The top row was plotted using RK4 and the bottom row RK45. It is clear that RK45 produces smoother plots and one explanation has to do with a larger discretization of the time interval. Whereas RK4 uses a for loop algorithm since the time step is constant for every iteration, RK45 uses a while loop because the algorithm readjusts the time step to minimize the error difference within a tolerance.

Above is a plot of comparison of computational time between RK4 and RK45 for solving Lorenz Attractor. Initial step size used is h=0.1s. Although, RK45 provides a more accurate estimate of the solution within a prescribed tolerance however it is more computationally expensive. Figure 6 shows the plot of comparison of the computational time between RK4 and RK45 for solving the Lorenz Attractor. Both methods start with the initial time step $h=0.1 seconds$. Horizontal axis is the time upper limit, t, in seconds; while the vertical axis is the runtime in seconds for each time upper limit. In light of the accuracy in data taking, for each time upper limit, the runtime was taken for three test runs and the average of these three runs are used to plot the relationship between time upper limit versus runtime. The graphs obtained from both methods are approximately linear line. By using the linear regression in MALAB, the slope of RK4 and RK45 are 0.00065 and 0.0035, respectively. This indicates that as time progresses the runtime for solving the Lorenz Attractor increases by a factor of 0.00065 and 0.0035 for using RK4 and RK45, respectively.

Lorenz Attractor with $\rho=28 \; \sigma=10 \; \beta=8/3$

The time spectrum on Figure 7 suggests that the system is always fluctuating and no equilibrium (dynamic or static) is ever reached. Moreover, for a deterministic system, one would expect Figure 6 to have a constant fluctuation (much like a sine wave or horizontal line), but it is clear that the attractor begins to exhibit a chaotic behavior after $t=13$. This is essentially the main idea of chaos theory - that deterministic system never achieves a steady state and trajectory never coincides with any other (See Figure 6). Using both Runge Kutta methods, one thing is certain, it doesn't matter how close two different initial conditions are, their trajectories will eventually diverge.

Conclusion

The fact that the Lorenz attractor is sensitive to its starting point to a plausible extent suggests that Weather is usually predictable only about a week ahead [11]. Long term prediction cannot be made due to the system’s chaotic behavior after some t time. A future study can be examined in the context of random data versus chaotic data. This page solely seeks a comparison between the two Runge Kutta methods as a numerical approximation to study the theory, but real life application calls for unwanted and corrupting noises. Theoretically, even when the system is deterministic, the real time series will consist of a few stochastic data, therefore one should be careful when defining the error as the difference between the time evolution of the ’test’ state and the time evolution of the nearby state. Finally, it is without a doubt that studying chaos theory will provide mankind better control of the predictibility of standard living, but we are far from comprehending its chaotic behavior. It is difficult to even a simplify a gaseous system in a small domain, let alone simulate the weather by breaking the atmosphere into many millions of interactive systems and predict the weather for next month since we don’t know the initial conditions at time t=0 very accurately.

References

1. Danforth, Christopher M. Chaos in an Atmosphere Hanging on a Wall, Mathematics of Planet Earth, 2013.
2. Fractal Foundation, What is Chaos Theory, Fractals are Smart: Science, Math and Art! http://fractalfoundation.org.
3. George, Yurkon T. (1997). Introduction to Chaos and It's Real World Applications. http://www.csuohio.edu/sciences/dept/physics/physicsweb/kaufman/yurkon/chaos.html.
4. Camhdag, Ezel. Runge Kutta Methods. Cankaya University, January 27, 2009.
5. Press, W. H. (1996). Chapter 16. Integration of Ordinary Differential Equations. Numerical recipes in fortran 77: the art of scientific computing (2nd ed., pp. 701-744). Cambridge [England: Cambridge University Press.
6. Press, William H. , Teukolsky, Saul A. , Vetterling, William T, and Flannery, Brian P. Numerical Recipes in C, Second Edition, Cambridge University Press 1992.
7. Bradley, Larry. Chaos and Fractals, 2010 http://www.stsci.edu/~lbradley/seminar/references.html.
8. Kellert, Stephen H. In the Wake of Chaos: Unpredictable Order in Dynamical Systems, University of Chicago Press 1993. p. 32.
9. Sparrow, Colin. The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors, Springer 1982.
10. Hilborn, Robert C. Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers (second ed.). Oxford University Press 2010.
11. Watts, Robert G. (2007). Global Warming and the Future of the Earth. Morgan and Claypool. p. 17.