Numerical solution of partial differential equations pdf

Since, from Eq. A mechanical device which applies a force propor- tional to velocity is a dashpot. The minus sign in this equation means that the force opposes the direction of motion, as required to be physically realizable. Thus, the application of this dashpot to the end of a finite length bar simulates exactly a bar of infinite length Fig.

This problem corresponds physically to two-dimensional steady-state heat conduction over a rectangular plate for which temperature is specified on the boundary. We attempt an approximate solution by introducing a uniform rectangular grid over the domain, and let the point i, j denote the point having the ith value of x and the jth value of y Fig.

Then, using central finite difference approximations to the second derivatives Fig. For example, consider the mesh shown in Fig. Although there are 20 mesh points, 14 are on the boundary, where the temperature is known. Thus, the resulting numerical problem has only six degrees of freedom unknown variables. The application of Eq. This linear system of six equations in six unknowns can be solved with standard equation solvers. Because the central difference operator is a 5-point operator, systems of equations of this type would have at most five nonzero terms in each equation, regardless of how large the mesh is.

Thus, for large meshes, the system of equations is sparsely populated, so that sparse matrix solution techniques would be applicable. Since the numbers assigned to the mesh points in Fig. Some equation solvers based on Gaussian elimination oper- ate more efficiently on systems of equations for which the nonzeros in the coefficient matrix are clustered near the main diagonal. Such a matrix system is called banded. Systems of this type can also be solved using an iterative procedure known as relaxation, which uses the following general algorithm: 1.

Initialize the boundary points to their prescribed values, and initialize the interior points to zero or some other convenient value e. Loop systematically through the interior mesh points, setting each interior point to the average of its four neighbors. Continue this process until the solution converges to the desired accuracy. For example, consider again the problem of the last section but with a Neumann, rather than Dirichlet, boundary condition on the right side Fig. The procedure is used in a variety of applications, including structural mechanics and dy- namics, acoustics, heat transfer, fluid flow, electric and magnetic fields, and electromagnetics.

Although the main theoretical bases for the finite element method are variational principles and the weighted residual method, it is useful to consider discrete systems first to gain some physical insight into some of the procedures. We let u2 and u3 denote the displacements from the equilibrium of the two masses m2 and m3.

The stiffnesses of the two springs are k1 and k2. This approach would be very tedious and error-prone for more complex systems involving many springs and masses. To develop instead a matrix approach, we first isolate one element, as shown in Fig. In the absence of the masses and constraints, this system is shown in Fig. This matrix corresponds to the unconstrained system.

Notice that, from Eqs. Without constraints, K is singular and the solution of the mechanical problem is not unique because of the presence of rigid body modes. K is symmetric. This property is a special case of the Betti reciprocal theorem in mechanics. An off-diagonal term is zero unless the two points are common to the same element. Thus, K is sparse in general and usually banded.

K is singular without enough constraints to eliminate rigid body motion. For spring systems, that have only one DOF at each point, the sum of any matrix column or row is zero. The forces must sum to zero, since the object is in static equilibrium.

We summarize the solution procedure for spring systems: 1. Generate the element stiffness matrices. Assemble the system K and F. Apply constraints. Compute reactions, spring forces, and stresses. However, matrix assembly for a truss structure a structure made of pin-jointed rods differs from that for a collection of springs, since the rod elements are not all colinear e.

Thus the elements must be transformed to a common coordinate system before the element matrices can be assembled into the system stiffness matrix. A typical rod element in 2-D is shown in Fig.

The four DOF for the element are also shown in Fig. These four values complete the first column of the matrix. Later, in the discussion of matrix transformations, we will derive a more convenient way to obtain this matrix. By using matrix partitioning, we can treat nonzero constraints and recover the forces of constraint. The grid point forces Fs at the constrained DOF consist of both the reactions of constraint and the applied loads, if any, at those DOF.

For example, consider the beam shown in Fig. When the applied load is distributed to the grid points, the loads at the two end grid points would include both reactions and a portion of the applied load. However, for general purpose software, such a capability is essential. An alternative approach to enforce this constraint is to attach a large spring k0 between DOF i and ground a fixed point and to apply a force Fi to DOF i equal to k0 U.

This spring should be many orders of magnitude larger than other typical values in the stiffness matrix. A large number placed on the matrix diagonal will have no adverse effects on the conditioning of the matrix. The large spring approach can be used for any system of equations for which one wants to enforce a constraint on a particular unknown. The main advantage of this approach is that computer coding is easier, since matrix sizes do not have to change.

Choose the large spring k0 to be, say, 10, times the largest diagonal entry in the unconstrained K matrix. For each DOF i which is to be constrained zero or not , add k0 to the diagonal entry Kii , and add k0 U to the corresponding right-hand side term Fi , where U is the desired constraint value for DOF i.

In two dimensions, we designate the four DOF associated with flexure as shown in Fig. Rotations are expressed in radians. For this element, note that there is no coupling between axial and transverse behavior. Transverse shear, which was ignored, can be added. The three-dimensional counterpart to this matrix would have six DOF at each grid point: three translations and three rotations ux , uy , uz , Rx , Ry , Rz.

In addition, for bending in two different planes, there would have to be two moments of inertia I1 and I2 , in addition to a torsional constant J and possibly a product of inertia I For most problems of interest in engineering, exact stiffness matrices cannot be derived. This figure also shows the domain modeled with several triangular elements.

A typical element is shown in Fig. Note that the number of undetermined coefficients equals the number of DOF 6 in the element. The linear approxima- tion in Eq. At the vertices, the displacements in Eq.

A is positive for counterclockwise ordering and negative for clockwise ordering. Note that, for this element, the strains are independent of position in the element. For this element, the stresses are constant independent of position in the element. We now derive the element stiffness matrix using the Principle of Virtual Work. Consider an element in static equilibrium with a set of applied loads F and displacements u. Note that, since the material property matrix D is symmetric, K is symmetric.

For example, in the finite element method, it is frequently more convenient to derive element matrices in a local element coordinate system and then transform those matrices to a global system Fig.

Transformations are also needed to transform from other orthogonal coordinate systems e. An orthonormal basis is defined as a basis whose basis vectors are mutually orthogonal unit vectors i. If we take the dot product of both sides of Eq. A matrix R satisfying Eq. That is, an orthogonal matrix is one whose inverse is equal to the transpose. R is sometimes called a rotation matrix. For example, for the coordinate rotation shown in Fig.

We recall that the determinant of a matrix product is equal to the product of the deter- minants. Also, the determinant of the transpose of a matrix is equal to the determinant of the matrix itself. That is, the square of the length of a vector is the same in both coordinate systems. A tensor of rank 0 is a scalar a quantity which rule v is unchanged by an orthogonal coordinate transformation. We now introduce tensors of rank 2. We now prove that cijkl is a tensor of rank 4.

We can write Eq. In this equation, u and F contain several subvectors, since u and F are the displacement and force vectors for all grid points, i. Consider the rod shown in Fig. This result agrees with that found earlier in Eq. For an isotropic material, cijkl must be an isotropic tensor of rank 4, thus implying at most three material constants on the basis of tensor analysis alone. In general, a functional is a function of a function.

More precisely, a functional is an integral which has a specific numerical value for each function that is substituted into it. The second term in Eq. This type of boundary condition is called a geometric or essential boundary condition. Consider the curve in Fig. Thus, the Euler-Lagrange equation, Eq. Assume that the bead starts from rest. From Eq. Also, from Eq. To solve Eq. The resulting cycloids for several end points are shown in Fig.

This brachistochrone solution is valid for any end point which the bead can reach. Thus, the end point must not be higher than the starting point. The end point may be at the same elevation since, without friction, there is no loss of energy during the motion. Thus, to enforce the constraint in Eq.

Note that, in 3-D, the integration is a volume integral. In general, a key conclusion of this discussion of variational calculus is that solving a differential equation is equivalent to minimizing some functional involving an integral.

It can be shown that the generalization of the Euler-Lagrange equation, Eq. For linear equa- tions, finite element solutions are often based on variational methods. The range is 1,2,3 in 3-D and 1,2 in 2-D. Let f x1 , x2 , x3 be a scalar function of the Cartesian coordinates x1 , x2 , x3. Using the comma notation and the summation convention, a variety of familiar quantities can be written in compact form. Given the functional, it is generally easy to see what partial differential equation corresponds to it.

The harder problem is to start with the equation and derive the variational principle i. To simplify this discussion, we will consider only scalar, rather than vector, problems. A scalar problem has one dependent variable and one partial differential equation, whereas a vector problem has several dependent variables coupled to each other through a system of partial differential equations. Let us subdivide the domain into triangular finite elements, as shown in Fig.

A typical element, shown in Fig. The number of DOF for an element represents the number of pieces of data needed to determine uniquely the solution for the element. The determinant 2A is positive or negative de- pending on whether the grid point ordering in the triangle is counter-clockwise or clockwise, respectively.

Since Eq. Notice that all the subscripts in this equation form a cyclic permutation of the numbers 1, 2, 3. That is, knowing N1 , we can write down N2 and N3 by a cyclic permutation of the subscripts. In general, the functions Ni are called shape functions or interpolation functions and are used extensively in finite element theory. In general, Ni is the shape function for Point i. We note that, since the shape functions for the 3-point triangle are linear in x and y, the gradients in the x or y directions are constant.

The displacement u x of an axial structural member a pin-jointed truss member Fig. For a linear variation of displacement along the length, it follows that Ni must be a linear function which is unity at Point i and zero at the other end.

We consider a slight generalization of the problem of Eq. The difference between this problem and that of Eq. The h term could arise in a variety of physical situations, including heat transfer due to convection where the heat flux is proportional to temperature and free surface flow problems where the free surface and radiation boundary conditions are both of this type.

Assume that the domain has been subdivided into a mesh of triangular finite elements similar to that shown in Fig. The functional which must be minimized for this boundary value problem is similar to that of Eq. We thus can consider a typical element. Thus, the last two terms in the functional make a nonzero contribution to the functional only if an element abuts S2 i. The two terms in F correspond to the body force and Neumann boundary condition, respectively.

The superscript e in the preceding equations indicates that we have computed the con- tribution for one element. We must combine the contributions for all elements. Consider two adjacent elements, as illustrated in Fig. In that figure, the circled numbers are the element labels. Because of the historical developments in structural mechanics, K is sometimes referred to as the stiffness matrix. Note also in Eq.

The zero gradient boundary condition is the natural boundary condition for this formulation, since a zero gradient results by default if the boundary is left free. To compute K from Eq. However, bi and ci are computed from Eq. Thus, two elements that are geometrically congruent and differ only by a translation will have the same element matrix. To compute the right-hand side contributions as given in Eq. For nonuniform f , the integral could be computed using natural or area coordinates for the triangle [7].

It is also of interest to calculate, for the linear triangle, the right-hand side contributions for the second term of Eq. For the second term of Eq.

To calculate H using Eq. Consider some triangular elements adjacent to a boundary, as shown in Fig. Since the integration in Eq. Thus, using the boundary shape functions of Eq. The first term for I is a quadratic form. In solid mechanics, I corresponds to the total potential energy. The first term is the strain energy, and the second term is the potential of the applied loads. Since strain energy is zero only for zero strain corresponding to rigid body motion , it follows that the stiffness matrix K is positive semi-definite.

For well-posed problems which allow no rigid body motion , K is positive definite. Three consequences of positive-definiteness are 1. All eigenvalues of K are positive. The matrix is non-singular. Gaussian elimination can be performed without pivoting. It was also assumed that the integral evaluated over some domain was equal to the sum of the integrals over the elements.

We wish to investigate briefly the conditions which must hold for this assumption to be valid. That is, what condition is necessary for the integral over a domain to be equal to the sum of the integrals over the elements?

Thus, we conclude that, for any functional of interest in finite element analysis, the integrand may be discontinuous, but we do not allow singularities in the integrand. This requirement is referred to as the compatibility requirement or the conforming requirement. For the 3-node triangular element already formulated, the shape function is linear, which implies that, in Fig. A nonconforming element would result in a discontinuous displacement at the element boundaries i. However, note that having displacement continuous implies that the displacement gradients which are proportional to the stresses are discontinuous at the element boundaries.

This property is one of the fundamental characteristics of an approxi- mate numerical solution. If all quantities of interest e. Thus, the approximation inherent in a displacement-based finite element method is that the displacements are continuous, and the stresses displacement gradients are discontinuous at element boundaries.

The best approximate solution will be one which in some sense minimizes R at all points of the domain. This is an example of a parabolic equation. We assume that the solution can be written as the product of a function of x and a function of t, ie.

Our discretised PDE equation 2. Solving equation 2. Figure 2. You can check that using the matlab code ForwardEuler. The trick with the Method of Lines is that it replaces all spatial derivatives with finite differences but leaves the time derivatives.

It is then possible to use a stiff ordinary differential equation solver on the time derivatives in the resulting system. So using equation 3. The matlab code for this example is BackwardEuler. Initial condition for Temperature distribution Variation of Temperature distribution with time using Backward Euler method 3 2. The solution using the Backward Euler method in figure 3. So even though there is more work in each time step inverting a matrix using the implicit Backward Euler method it allows larger time steps than the explicit Forward Euler method.

How can we speed up this calculation? By using iterative methods. However they are ideal for finite difference methods involving solution of large sparse matrices. Different methods have different convergence times and for big inverse matrix problems are much faster than direct matrix inverse methods.

The Jacobi method needs 0 17 iterations to converge so takes nearly twice as long as the Gauss-Seidel method. Because A should be used to solve for U is sparse and diagonally dominant, iterative solutions are ideal here. However we see in the next section 5. These are solved using LU decomposition. For more details see Schilling and Harris, p.

Because the problem has circular symmetry ie. We solve equation 5. Initial condition for Temperature distribution Temp after 1 year If we let:! We assume that independent solutions eigenmodes of equation 7. Equation 7. We will show this in the next section.

The matlab code is Lax Flux. The numerical solution changes for different time steps depending on whether or not the scheme is stable or dispersion is present. We plot the numerical solution in figure 8.

Figure 8. We consider friction due to viscosity of medium and density of string. Again we use central difference for Uxx and Utt as in section 8. The matlab code is Wave1DFriction. Figure 9. A weak solution of 9.

Then 9. To solve equation 9. We solve 9. Using FEM we seek an approximate solution to 9. We will show that the stiffness matrix C is tridiagonal for this example. The weak solution for U satisfies the variational form for Equation 9. Similarly to the 1-D case we seek an approximate solution to equation 9. The 2D hat function is plotted in figure 9. We require that equation 9. Equation 9. However we still need to approximate time derivatives.