Search in book...
Toggle Font Controls
Create new playlist

Name your new playlist

Playlist description (optional)
Sign In

Email address

Password

Forgot Password?

or

Continue with Facebook

Continue with Google
Sign Up

Full Name

Email address

Confirm Email Address

Password

or

Continue with Facebook

Continue with Google

4
Symbolic Operations

4.1 Introduction

One part of solving engineering problems is to apply one or more solution methods until a form that is amenable to numerical procedures and computer evaluation is obtained. We shall illustrate in this chapter how one can use the symbolic capabilities of Mathematica to meet these goals. One advantage of using symbolic solutions is that once their symbolic solution is known it can be converted to a function and used to perform parametric studies or used in other applications.

We shall illustrate two types of symbolic manipulations. The first type is concerned with the simplification and manipulation of the symbolic expressions to attain a form more suitable for one’s end usage. The second type is to perform a mathematical operation on a symbolic expression such as finding its derivative, integrating it, and the like.

Several simplification and manipulation commands that are typically used on symbolic expressions are listed in Table 4.1. Also included in this table are illustrations of the usage of these functions.

images — **Table 4.1** Manipulation of symbolic expressions

The symbolic mathematical operations that will be illustrated are:

Solving equations—Solve[]
Limits—Limit[]
Power series—Series[] and Coefficient[]
Optimization—Maximize[]/Minimize[]
Differentiation—D[]
Integration—Integrate[]
Solutions to ordinary differential equations—DSolve[]
Solutions to some partial differential equations—DSolve[]
Laplace transform—LaplaceTransform[] and InverseLaplaceTransform[]

Not all differential equations and integrals have symbolic solutions; in these cases, one must obtain their solutions numerically. Therefore, in the next chapter, we shall re-visit some of these same operations using functions that employ numerical procedures to obtain their solutions.

For many uses of the symbolic operations, the most compact final forms of the results are often attained by operating on individual collections of terms and then placing each of them in their most compact form. Thus, the reader is encouraged to interact with the symbolic capabilities of Mathematica one line of code at a time. Several of the examples that follow employ this strategy.

4.2 Assumption Options

Many of the functions that we shall discuss in this chapter, including many of those appearing in Table 4.1, allow one to include assumptions concerning one or more of the symbols appearing in the expression being evaluated or manipulated. In those functions that accept Assumptions, the form can be either


…[…,Assumptions->{asum1, asum2,…}]


…[…,Assumptions->asum1 && asum2 && …]

where asumN are any of the assumptions appearing in the right-hand column of Table 4.2.

Table 4.2 Common assumptions that can be used in those Mathematica functions that accept them

Assumption*	Mathematica expression for assumption:⁺ …,Assumptions->…&&… or …,Assumptions->{…,…}
x is a real number	Element[x,Reals] or Reals
x is an integer	Element[x,Integers] or Integers
x is a complex number	Element[x,Complexes] or Complexes
x is a prime number	Element[x,Primes] or Primes
x > n or x ≥ n	x>n or x>=n
x < n or x ≤ n	x<n or x<=n
m < x < n or m ≤ x < n or m < x ≤ n or m ≤ x ≤ n	m<x<n or m<=x<n or m<x<=n or m<=x<=n

^*x = variable name or a list of variable names that appear in the expression being operated on by the function; n, m = numerical or symbolic values.

⁺The symbol can be found in the Special Characters palette by selecting the Symbols tab, then the ≠ tab, and then going to the right-most column of the next to last row.

4.3 Solutions of Equations: Solve[]

The symbolic solutions to a set of equations, if they exist, are obtained by using either


sol=Solve[{eq1L==eq1R,eq2L==eq2R,…},vars]


sol=Solve[eq1L==eq1R && eq2L==eq2R &&…,vars]

where eqNR can be zero and vars is a list of the solution variables. The number of solution variables equals the number of equations.

Example 4.1

Solution of Nonlinear Equations

Consider the two equations

Their solution is obtained by using


xysol=Solve[x^2+y^2==1 && x+y==a,{x,y}]

which yields

It is seen that xysol is a 2×2 array.

We shall now detail the procedure that can be used to access the two sets of solutions. Then, a compact syntax will be introduced that we shall use frequently in this chapter and in subsequent chapters. A typical way that these results are accessed is to place them in a list. This is done with


z={x,y}/.xysol

which displays

Therefore, one solution is (x₁, y₁) = (z[[1,1]],z[[1,2]]) and the other solution is (x₂, y₂) = (z[[2,1]],z[[2,2]]).

A more compact way to get these results into an accessible list form is by combining the above syntax as follows


z={x,y}/.Solve[x^2+y^2==1 && x+y==a,{x,y}]

To determine the values of (x₁, y₁) and (x₂, y₂) for, say, a = 0.5, we use the substitution rule as follows


v=z/.a->0.5

which displays


{{−0.411438,0.911438},{0.911438,−0.411438}}

Thus, for a = 0.5, x₁ = v[[1,1]] = −0.411438, and y₁ = v[[1,2]] = 0.911438 and x₂ = v[[2,1]] = 0.911438, and y₂ = v[[2,2]] = −0.411438.

Example 4.2

Solution of a System of Equations

In Section 2.5, a solution to the following system of equations was obtained by using LinearSolve

(a)

We shall change this system of equations slightly by replacing the value on the right-hand side of the second equation to an unknown constant b. Then, Eq. (a) becomes

(b)

Using the compact syntax of the previous example, the solution to Eq. (b) is obtained with


x₁₂₃={x₁,x₂,x₃}/.Solve[2 x₁+3 x₂+5 x₃==−1 &&
7 x₁+11 x₂+13 x₃==b && 17 x₁+19 x₂+23 x₃==31,{x₁,x₂,x₃}]

where the appropriate template from the Basic Math Assistant palette was used. The execution of this statement yields

To place these results in a slightly more convenient form, we use


x₁₂₃=Flatten[x₁₂₃]

which gives

If an annotated display of the results were desired, the following is used


 Print["x₁=",x₁₂₃ [[1]]]
 Print["x₂=",x₁₂₃ [[2]]]
 Print["x₃=",x₁₂₃ [[3]]]

which displays

To show that these values agree with those obtained in Section 2.5, we set b = 7.5. Thus,


x₁₂₃= x₁₂₃/.b->7.5

yields the values obtained in Section 2.5; that is,


{3.9359,0.942308,−2.33974}

We shall again obtain the solution of Eq. (a), but this time we shall assume that the numbers have units. In particular, it is assumed that the values multiplying x_n have the units of N·m⁻¹ and those on the right-hand side of the equal sign have the units of N. In this case, we use Section 1.11, the compact syntax introduced in Example 4.1, and the above results to obtain


nm="Newtons"/"Meters"; ns="Newtons";
 x₁₂₃=Flatten[{x₁,x₂,x₃}/.Solve[
  {Quantity[2,nm] x₁+Quantity[3,nm] x₂+
   Quantity[5,nm] x₃==Quantity[−1,ns]&&
  Quantity[7,nm] x₁+Quantity[11,nm] x₂+
   Quantity[13,nm] x₃==Quantity[7.5,ns]&&
  Quantity[17,nm] x₁+Quantity[19,nm] x₂+
   Quantity[23,nm] x₃==Quantity[31,ns]},{x₁,x₂,x₃}]]/.b->7.5;
 Print["x₁=",x₁₂₃[[1]]]
 Print["x₂=",x₁₂₃[[2]]]
 Print["x₃=",x₁₂₃[[3]]]

which displays


 x₁=3.935 m
 x₂=0.942308 m
 x₃=−2.33974 m

4.4 Limits: Limit[]

The limit of an expression is obtained from


Limit[expr,x->xo,Assumptions->…]

where expr is a function of x and xo is a symbolic quantity or a numerical quantity, where -Infinity < xo < Infinity.

For example, to determine

we use


Limit[(1−Cos[x])/x,x->0]

and obtain 0.

Example 4.3

Heaviside Theta Function

The Heaviside theta function H(x) is defined as follows: H(x) = 1, x > 0 and H(x) = 0, x < 0. It is undefined for x = 0. One mathematical construct that has the properties of the Heaviside theta function is (see also Example 4.20)

To verify this, we determine the limits for ± x, x ≠ 0. Thus,


hxp=Limit[1/2+1/π ArcTan[x/t],t->0]
hxm=Limit[1/2+1/π ArcTan[−x/t],t->0]

which gives

To simplify these expressions, we use PowerExpand. Thus,


 hxp=PowerExpand[hxp]
 hxm=PowerExpand[hxm]

which gives that hxp = 1 and hxm = 0.

Example 4.4

Limit Using Assumptions

As indicated above, the Limit function also accepts Assumptions. To see what difference this can make on the final form of the limit, consider the limit

If we make no assumptions about n, then the limit is obtained from


Limit[Sin[n x]/Sin[2 n x],x->π]

which displays

On the other hand, if n is an integer we can rewrite the limit statement as


Limit[Sin[n x]/Sin[2 n x],x->π,Assumptions->n  Integers]

which gives

This is the same result that one obtains from the previous result, except an additional step is necessary.

Example 4.5

Natural Frequencies of Beams

The limit process can also be used to obtain special cases of very general solutions to engineering applications. Consider the following characteristic equation in terms of a frequency parameter Ω for an Euler beam that is supported at each end by translational springs expressed by the nondimensional constants 0 ≤ K_L ≤ ∞ and 0 ≤ K_R ≤ ∞ and torsional springs expressed by the nondimensional constants 0 ≤ T_L ≤ ∞ and 0 ≤ T_R ≤ ∞

where

and

Special cases of this general result can be obtained by selecting the appropriate limiting values for the spring constants. Thus, for the case of a beam clamped at both ends, K_L → ∞, K_R → ∞, T_L → ∞, and T_R → ∞. This is equivalent to letting a_1n → ∞, a_2n → ∞, b_1n → ∞, and b_2n → ∞. Thus, we divide the characteristic equation sequentially by these quantities and after each division we take the limit to arrive at the characteristic equation for a beam clamped at both ends.

The program is


 z₁=b₁ b₂ (a₁+a₂)+b₁-b₂;
 z₂=a₁ a₂ (b₁-b₂)-a₁-a₂;
 z₃=a₁ a₂+b₁ b₂;
 z₄=1-a₁ a₂ b₁ b₂;
 z₅=a₂ b₂-a₁ b₁;
 z₆=a₁ b₂-a₂ b₁;
 ce=z₁ (Cos[Ω] Sinh[Ω]+Sin[Ω] Cosh[Ω])+2 z₆ Cos[Ω] Cosh[Ω]+
  z₂ (Cos[Ω] Sinh[Ω]−Sin[Ω] Cosh[Ω])−2 z₃ Sin[Ω] Sinh[Ω]+
  z₄ (Cos[Ω] Cosh[Ω]−1)+z₅ (Cos[Ω] Cosh[Ω]+1);
 p=Limit[Limit[Limit[Limit[ce/a₁,a₁−>∞]/a₂,a₂−>∞]/b₁,b₁−>∞]/
  b₂,b₂−>∞]

which, upon execution, gives


1−Cos[Ω] Cosh[Ω]

Example 4.6

Determination of the Poles of an Expression

In using the inversion theorem for Laplace transforms, it is necessary to know whether a pole (zero of a denominator) is a simple pole; that is, the expression for the pole of the denominator appears as (z − z_o), where z_o is the pole. If the denominator appears as (z − z_o)^m it is called an mth-order pole.

Consider the expression

To determine if z = 0 is a pole, we take the limit of f(z) as z → 0. Thus,


Limit[Exp[z t] (η z Cosh[z]−Sinh[η z])/(z^3 Cosh[z]),z−>0]

gives

Thus, z = 0 is not a pole; if it were a pole, the limit would have been ∞.

4.5 Power Series: Series[], Coefficient[], and CoefficientList[]

One can obtain an nth-order power series approximation of an expression about the point x = x_o with


Series[f,{x,xo,n}]

where f is the expression for which the series is to be obtained.

To illustrate the use of Series, we shall get a fourth-order series of the cosine of x about the point a. The statement to do this is


p=Series[Cos[x],{x,a,4}]

which gives

The expression O[x-a]⁵ represents a term on the order of (x − a)⁵. When a series is given in this form, one cannot use it for further numerical evaluation; that is, for example, p/.x−>2 results in an error message. To obtain an expression without this term and in a form that permits further numerical evaluation, the function Normal is used as follows


p=Normal[Series[Cos[x],{x,a,4}]]

Execution of this command results in the previous result, but without the last term O[x-a]⁵.

To be able to compare this result to those that follow, we shall expand p as follows


Expand[p]

which results in

The coefficients of the series can be accessed two ways. If expans is a polynomial in the powers of (x − a) up to the nth power, then the first way that the coefficients can be accessed is by using


Coefficient[expans,(x-a),m]]

which retrieves the coefficient associated with (x − a)^m, n ≥ m ≥ 1, or with


CoefficientList[expans,x]

which gives a list of all the coefficients up to and including xⁿ⁻¹ (n ≥ 1).

For example, continuing with the expansion of cosine of x above,


Coefficient[p,(x-a),2]]

gives

whereas


Coefficient[p,x,2]]

gives

On the other hand,


CoefficientList[p,x]

displays

This list contains the coefficients of the expansion shown in Eq. (4.1), where the first element in the list contains the terms multiplying x⁰, the second element the terms multiplying x¹, and so on.

Example 4.7

Perturbation Solution #1

The solution to

with the initial conditions of y(0) = 0 and dy(0)/dt = 1 is

It is assumed that 1 and the expansion of this solution up to terms including ² is sought. The statement to do this is


Normal[Series[(1+) (1-Exp[− t])/^2−t/,{,0,2}]]

which displays

Example 4.8

Perturbation Solution #2

Perturbation theory can be used to obtain the approximate roots of transcendental equations. For example, consider the expression

(a)

where μ is a small quantity. To determine the value of x that satisfies this equation using the perturbation technique, one assumes a solution of the form

(b)

where X_n is to be determined.

If we assume that N = 3 and substitute Eq. (b) into Eq. (a), we have

To obtain the values of X_n, the above expression is expanded by expressing the exponential function as a series and then by collecting the terms multiplying like powers of μ. These terms are then set equal to zero and a system of equations is obtained from which X_n can then be determined.

We shall introduce the solution process in stages. First, we determine the coefficients of the powers of μ using


trans=Normal[(X₀+μ X₁+μ^2 X₂+μ^3 X₃)^2−1
  −μ Series[Exp[X₀+μ X₁+μ^2 X₂+μ^3 X₃],{μ,0,2}]]

which gives

Next, we create a list of equations in the format required by Solve. First, we create a list of the coefficients of μ as follows


zz=CoefficientList[trans,μ]

This gives

We now convert the elements of this list into equations equaling zero by appending “==0” to each of these coefficients. Thus,


q=Map[#==0 &,zz]

which produces

Lastly, we solve these four equations appearing in q for X₀ to X₃ using Solve as follows


Solve[q,{X₀,X₁,X₂,X₃}]

which gives

Then using the above result and Eq. (b), the two roots are approximated by

4.6 Optimization: Maximize[]/Minimize[]

For a function f(x), the value of x = x_m at which this function is a maximum/minimum and the magnitude of f(x_m) is obtained from


 Maximize[{f,con},x]
 Minimize[{f,con},x]

where f is the function to be maximized or minimized, con are any constraints that f is subject to, and x is the independent variable. The output of these functions is often fairly complex so that several other functions are frequently used in conjunction with them.

To isolate the expression corresponding to the value of f(x_m), Maximize and Minimize, respectively, are replaced with


 MaxValue[{f,con},x]
 MinValue[{f,con},x]

and to isolate the expression corresponding to x_m, Maximize and Minimize, respectively, are replaced with


 ArgMax[{f,con},x]
 ArgMin[{f,con},x]

One additional function that is used to further reduce the results to a specific region is


Refine[exp,assum]

where exp is a symbolic expression, in this case the output from one of the above functions, and assum are the assumptions that specify the region of interest.

The use of these functions is illustrated in the following example.

Example 4.9

Peak Amplitude Response of a Single Degree-of-Freedom System

The amplitude response function for a single degree-of-freedom system is given by

where Ω is a frequency ratio and 0 < ζ < 1. Expressions for the maximum value and the value of Ω at which this maximum occurs are obtained from


homega[Ω_,ζ_]:=1/Sqrt[(1−Ω^2)^2+(2 ζ Ω)^2]
 hmax=Maximize[{homega[Ω,ζ],0<ζ<1,Ω>0},Ω]//Simplify

which displays

It is seen that the solution of interest is that given for the region 0 < ζ < . Hence, the elements of this solution can be accessed in the following manner. To isolate the expression for H(Ω_m), we instead obtain the solution with


hmax=MaxValue[{homega[Ω,ζ],0<ζ<1,Ω>0},Ω]//Simplify;
 sh=Refine[hmax,0<ζ<1/Sqrt[2]]

which displays

Further simplification is done by using


hmax=PowerExpand[Factor[sh],Assumptions−>0<ζ<1/Sqrt[2]]

which displays

To obtain Ω_m, we re-solve the expression using


omegmax=ArgMax[{homega[Ω,ζ],0<ζ<1,Ω>0},Ω]//Simplify;
 so=Refine[omegmax,0<ζ<1/Sqrt[2]][[1]]
 Solve[so[Ω]==0,Ω]

which results in

Thus,

and, therefore,

4.7 Differentiation: D[]

Differentiation with respect to one or more variables is performed using


D[f,{x1,n1},{x2,n2},…]

where f is an expression in one or more of the variables {x1,x2,…} and n1,n2,… are the orders of the derivatives in each variable. Thus, D also obtains the partial derivative. There are several equivalent syntaxes for D that, along with examples of their usage, are shown in Table 4.3. Also presented in this table are examples of D applied to several classes of functions.

For an example, consider the expression

To determine we proceed as follows


D[x^3 y^4,{x,2},{y,3}] (* or D[x^3 y^4,x,x,y,y,y] *)

which gives


144 x y

When f is a function of one variable, D can be replaced by a prime (′) provided that f is explicitly represented as a function of the independent variable. Consider the following expression

Then the second derivative of y can be determined using primes as follows


y[x_]=x^3 Cos[a x];
y’’[x] (* or D[y[x],x,x] *)

The execution of these statements displays


6 x Cos[a x]-a² x³ Cos[a x]−6 a x² Sin[a x]

The command D can also be used to perform differentiation using functional representation and a change of variables. Consider the function f(u), where u = z(x). Then df/dx is obtained from


u=z[x];
 D[f[u],x] (* or in a compact form as D[f[z[x]],x] *)

which displays


f’[z[x]] z’[x]

where the prime denotes the derivative with respect to its argument; that is,

For example, if z = bx^g, then we can use the previous command and the transformation rule to determine df/dx as follows


D[f[z[x]],x]/.D[z[x],x]->D[b x^g,x]

which gives


b g x^−1+g f’[z[x]]

An expression for d²f/dx² is obtained from


D[f[z[x]],x,x]

which displays


z’[x]^2 f’’[z[x]]+f’[z[x]] z’’[x]

where the prime denotes the derivative with respect to its argument; that is,

Again letting z = bx^g, we can use the previous command and the transformation rule to determine d²f/dx² as follows


Simplify[D[f[z[x]],x,x]/.{D[z[x],x]->D[b x^g,x],
D[z[x],x,x]->D[b x^g,x,x]}]

which yields

We shall now illustrate the use of D with several additional examples.

Example 4.10

Radius of Curvature

The curvature of a plane curve in terms of the parametric quantities x(t) and y(t) is given by

where the prime indicates the derivative with respect to t. We shall determine κ when

where a is a constant.

The program to determine κ is


xp=D[a (t-Tanh[t]),t];
yp=D[a Sech[t],t];
kap=Simplify[(xp D[yp,t]-yp D[xp,t])/(xp^2+yp^2)^(3/2)]

Execution of these instructions gives

This result can be further simplified by using PowerExpand. Thus, we add the statement


kap=PowerExpand[kap]

which, upon execution, gives the final result

Example 4.11

Euler–Lagrange Equation

The governing differential for the dynamic response of a thin beam resting on an elastic foundation can be obtained from the Euler–Lagrange equation

where F is an expression relating the kinetic energy, the potential energy, and the external work of the system and

The quantity F for this beam system is given by

where c_k are constants.

We shall now use these equations to derive the governing equation as follows


 fF=1/2 c₁ D[w[x,t],t]^2−1/2 c₂ D[w[x,t],{x,2}]^2
  −1/2 c₃ w[x,t]^2+g[x,t] w[x,t];
 eqn=D[fF,w[x,t]]+D[D[fF,D[w[x,t],x,x]],x,x]
  −D[D[fF,D[w[x,t],t]],t]

Upon execution, we obtain


g[x, t]min usc₃ w[x, t]min usc₁ w^{(0, 2)}[x, t]min usc₂ w^{(4, 0)}[x, t]$? >

which is the governing equation of the beam. The superscripts in this expression are interpreted as follows. The first value in the superscript indicates the order of the derivative of the first variable in w[x,t], x in this case, and the second superscript indicates the order of the derivative of the second variable in w[x,t], t in this case. Thus, in standard form, this equation is, after rearrangement,

Variational Methods Package: The results in this example could also have been obtained with the Variational Methods Package. First, the package is loaded into the Mathematica session with


Needs["VariationalMethods'"]

Then, the governing equation is obtained with


EulerEquations[fF,w[x,t],{x,t}]

which yields the previous result.

4.8 Integration: Integrate[]

The indefinite integral with respect to one or more variables is given by


Integrate[expr,x,y,…,Assumptions->…]

where expr is, in general, a function of x, y, ….

The definite integral with respect to one or more variables is obtained using


 Integrate[expr,{x,xlow,xup},{y,ylow,yup},…,
   Assumptions->…]]

where expr is, in general, a function of x, y, …, xlow is the lower x limit and xup is the upper x limit and similarly for the y variable ylow is the lower y limit and yup is the upper y limit. Either or both of the limits can be infinity; that is, Infinity (or ∞ from the appropriate template.)

When the integral cannot be found in terms of known functions, the system returns the Integrate statement as entered.

Before proceeding with several examples of the usage of Integrate, a command that has utility in evaluating integrals and in finding Laplace transforms, as well as other applications, is Piecewise, which is defined as


Piecewise[{{expr1,cond1},{expr2,cond2},…}]

where exprN are expressions and condN are most typically the ranges over which each exprN is to apply and is most often of the form a_N ≤ x ≤ b_N. An illustration of its usage is shown in Example 4.17.

Integration can also be performed using the appropriate template from the Advanced portion of the Calculator section of the Basic Math Assistant. Its usage is illustrated in Example 4.12.

Example 4.12

Integration #1

We shall evaluate the definite integral

where 0 ≤ ζ < 1. Then,


 Simplify[Exp[−ζ t]/Sqrt[1−ζ^2]*
   Integrate[Exp[ζ η] Sin[Sqrt[1−ζ^2] (t−η)],{η,0,t}]]

yields

Another way to program this function is to use the appropriate templates from the Basic Math Assistant. In this case, the above program becomes

which yields the previously obtained result.

Example 4.13

Fourier Coefficients

We shall determine the value of

where

and k, n, and m are integers.

We start by computing c_n for the case where n ≠ m. Then,


 an=Simplify[Integrate[Sin[m t] Cos[n t]/π,{t,0,2 π k/m}],
  Assumptions->{k,n}Integers];
 bn=Simplify[Integrate[Sin[m t] Sin[n t]/π,{t,0,2 π k/m}],
  Assumptions->{k,n}Integers];
 cn=PowerExpand[Simplify[Sqrt[an^2+bn^2]]]

which yields

For the case of n = m, we have


 an=Simplify[Integrate[Sin[n t] Cos[n t]/π,{t,0,2 π k/n}],
  Assumptions->{k,n}Integers];
 bn=Simplify[Integrate[Sin[n t]^2/π,{t,0,2 π k/n}],
  Assumptions->{k,n}Integers];
 cn=PowerExpand[Simplify[Sqrt[an^2+bn^2]]]

and obtain

Thus,

Example 4.14

Integration #2

Consider the following integral, which has infinity as one of its limits

Its evaluation is obtained with


Integrate[Exp[-a t^2],{t,0,∞},Assumptions->a>0]

which produces

Example 4.15

Integration #3

Integrals of special functions can also be determined. Consider the following integral that involves the Bessel function of the first kind of order m, m ≥ 0,

Its evaluation is obtained with


Integrate[r BesselJ[m,r/a]^2,{r,0,a},Assumptions->m>=0]

which gives

If the traditional form of the result is desired, then


 TraditionalForm[Integrate[r BesselJ[m,r/a]^2,{r,0,a},
    Assumptions->m>=0]

is used and the following is displayed

Example 4.16

Integration #4

As part of performing an inversion of a Laplace transform, the convolution integral leads to the following integral

The evaluation of this integral is obtained with


Integrate[Cos[η^2/(4 μ)]/(Sqrt[μ] Sqrt[(τ-μ)]),{μ,0,τ},
  Assumptions->η>0]/π

The execution of this statement displays

where FresnelC and FresnelS are the Fresnel integrals given by

Thus,

Example 4.17

Integration #5

Consider the integral

where

The evaluation of this integral is obtained with the aid of Piecewise as follows


 f=Piecewise[{{a₁,0<=x<b},{0,b<=x<d},{-a₂,d<=x<=c}}];
 Simplify[Integrate[f,{x,0,c},Assumptions->0<b<d<c]]

which gives


a₁ b+a₂ (−c+d)

Example 4.18

Cauchy Integral Formula

The Cauchy integral formula is

where g(z) has z_k, k = 1, 2, …, K, singular points all contained within C and Res stands for the residue of g(z) evaluated at each z_k. The singular points are those values of z_k for which g(z) → ∞. The residues are determined by using for each z_k


Residue[gz,{z,zk}]

where gz = g(z) and zk = z_k.

As an example, let

The values of z_k are determined from


zk=z/.Solve[z (z^2+9)==0,z]

which yields


{0,−3 I,3 I}

Then, the value of the integral is obtained from


 2 n I Total[Table[Residue[Cos[z]/(z (z^2+9)),{z,zk[[n]]}],
  {n,1,3}]]//Simplify

which displays

4.9 Solutions of Ordinary Differential Equations: DSolve[]

The symbolic solutions to both ordinary and partial differential equations are obtained with


DSolve[{eqn1,eqn2,…},{y1,y2,…},{x1,x2,…}]

where {eqn1,eqn2,…} is a list of coupled differential equations and their boundary conditions, {y1,y2,…} is a list of the dependent variables appearing in {eqn1,eqn2,…}, and {x1,x2,…} is a list of the independent variables. When the system of equations is composed of ordinary differential equations, there is only one independent variable. When the boundary conditions are omitted, the system obtains a solution in terms of unknown constants C[n], where n goes from 1 to the order of the differential equation. In the list {eqn1,eqn2,…}, the representation of each dependent variable has to be given explicitly as y1[x1,x2,…], y2[x1,x1,…], …. Not all equations will have a symbolic solution. However, the solutions to many classical equations and some variations thereof are available and these solutions are presented in terms of special functions where appropriate. Several examples are shown in Table 4.4.

We shall illustrate the use of DSolve with the following examples. In several of these cases, the results will be plotted using the basic form of Plot, which is given in Table 6.1.

Example 4.19

Solution to an Inhomogeneous Ordinary Differential Equation

Consider the following equation

where γ > 0 and the initial conditions are y(0) = y′(0) = 0.

The solution is obtained as follows


v=DSolve[{y’’[t]+y[t]==a Exp[-γ t],y[0]==0,y’[0]==0},y,t]

which yields

To access the solution and simplify it, we use


p=FullSimplify[y[t]/.v[[1]]]

which gives

This result can also be obtained by using the more compact syntax


 p=FullSimplify[y[t]/.DSolve[{y’’[t]+y[t]==a Exp[-t γ],
  y[0]==0,y’[0]==0},y,t][[1]]]

Additionally, if both y(t) and y′(t) were of interest, then the previous statement would be written as


 p=FullSimplify[{y[t],y’[t]}/.DSolve[{y’’[t]+y[t]==
   a Exp[-t γ],y[0]==0,y’[0]==0},y,t][[1]]]

which results in the two-element array

The first element p[[1]] = y(t) contains the result previously obtained and the second element p[[2]] = y′(t) contains the expression for the first derivative.

These results can be evaluated for plotting or evaluated at specific values of a, γ, and t. For example, if we are interested in y(t) evaluated at a = 1, t = 0.4, and γ = 1, 3, 5, 7, then because a, γ, and t are global variables


TableForm[Table[{γ,p[[1]]/.{a->1.,t->0.4}},{γ,1,7,2}],
  TableHeadings->{None,{"γ","y(0.4)"}}]

yields

The symbolic solution of the differential equation for, say, y(t) can be made into a function for subsequent use as follows


 uu[t_,γ_,a_]=p[[1]];

Then,


 uu[0.4,7,1]

displays


 0.0373135

which agrees with one of the values previously obtained. Notice that in the definition of uu there is no colon before the equal sign.

The expressions for y(t) and y′(t) can be plotted for a = 1.0 and γ = 1.0 as follows


 a=1.; γ=1.;
 Plot[p,{t,0,20},AxesLabel->{"t","y(t)"}]

which produces Figure 4.1. In this figure, the curve with the largest initial magnitude is y(t).

It is mentioned that p is only available during this notebook session. Therefore, if one wanted to use this result in another application where a and γ varied over a range of values and it was not computational efficient to solve the differential equation symbolically each time, then one can copy the symbolic result into another procedure and go from there. For example, using the previous result for the solution to the differential equation, one could create the following function in the current notebook or in another notebook

Figure 4.1 Graph of the numerical evaluation of the symbolic results of Example 4.19

Example 4.20

Analysis of Beams

The nondimensional equation of a thin beam subjected to a concentrated load of unit magnitude at x_o is given by

where δ(x) is the delta function. If it is assumed that the beam is hinged at each end, then the boundary conditions are

Noting that the delta function is implemented by using DiracDelta, the symbolic solution is obtained from


 hhbeam=Simplify[y[x]/.DSolve[{y’’’’[x]==DiracDelta[x-xo],
  y[0]==0,y’’[0]==0,y’’[1]==0,y[1]==0},y,x][[1]],
  Assumptions->1>xo>0]

The execution of this statement gives

where HeavisideTheta[x] is defined as 0 for x < 0 and 1 for x > 0. It is undefined at x = 0. This function differs from UnitStep[x], which is defined as 0 for x < 0 and 1 for x ≥ 0. In addition, the derivative of HeavisideTheta[x] ; that is, D[HeavisideTheta[x],x], gives DiracDelta. The derivative of UnitStep[x] does not reduce to DiracDelta.

We shall obtain a graph of y, y′, y′′, and y′′′ when xo = 0.7. Since we are interested in y and its derivatives, the above statement is rewritten as follows (recall Example 4.19)


 hhbeam=Simplify[{y[x],y’[x],y’’[x],y’’’[x]}/.DSolve[
   {y’’’’[x]==DiracDelta[x-xo],y[0]==0,y’’[0]==0,y’’[1]==0,
   y[1]==0},y,x][[1]],Assumptions->1>xo>0]

Hence, y(x) = hhbeam[[1]], y′(x) = hhbeam[[2]], and so on. We use GraphicsGrid to display the results on a 2×2 grid. See Table 6.17 for the use of GraphicsGrid. Then


 xo=0.7;
 plt[n_]:=Plot[hhbeam[[n]],{x,0,1}]
 GraphicsGrid[{{plt[1],plt[2]},{plt[3],plt[4]}}]

gives the results shown in Figure 4.2. The function plt was created to improve readability of the program.

Figure 4.2 Graph of the numerical evaluation of the symbolic results of Example 4.20

Example 4.21

Deformation of a Timoshenko Beam

Consider the following coupled second-order equations in terms of nondimensional quantities

which are valid over the region 0 ≤ x ≤ 1. These equations describe the static deformation of a Timoshenko beam. We shall assume that the boundary conditions are W(0) = W(1) = 0 and Ψ(0) = Ψ(1) = 0. Then, if we let q(x) = 1, the solutions are obtained from


 sol=Flatten[Simplify[{w[x],s[x]}/.DSolve[{w’’[x]-s’[x]==1,
   α s’’[x]-s[x]+w’[x]==0,w[0]==0,w[1]==0,s[0]==0,s[1]==0},
   {w[x],s[x]},x]]]

where we have employed the compact notation of Examples 4.19 and 4.20. The execution of this instruction gives

Therefore, the individual solutions are accessed with W(x) = sol[[1]] and Ψ(x) = sol[[2]]. If the value of W(0.3) when α = 2 is desired, then its value is obtained with


sol[[1]]/.{α->2.,x->0.3}

which yields −0.105919.

If a plot of the W(x) and Ψ(x) is desired for, say, α = 0.1, then the following is used


 α=0.1;
 Plot[sol,{x,0,1}]

Example 4.22

Logistic Equation

Consider the following nonlinear equation, called the logistic equation,

which has the initial condition x(0) = h. The solution is obtained from


 ss=x[t]/.DSolve[{x’[t]==x[t] (1-x[t]),x[0]==h},x[t],t][[1]]

which yields

To plot these results for h = 0.01, we use


 h=0.01;
  Plot[ss,{t,0,10},AxesLabel->{"t","x"}]

which generates Figure 4.3.

Figure 4.3 Graph of the numerical evaluation of the symbolic result of Example 4.22

Example 4.23

System of First-Order Equations and the Matrix Exponential

Consider the following system of first-order equations with constants coefficients

subject to the following initial conditions: x₁(0) = −2, x₂(0) = 1, and x₃(0) = 4. The solution is obtained with DSolve as follows.


 Expand[Flatten[{x₁[t],x₂[t],x₃[t]}/.DSolve[
  {x₁’[t]==6 x₁[t]+3 x₂[t]−2 x₃[t],
  x₂’[t]==−4 x₁[t]- x₂[t]+2 x₃[t],
  x₃’[t]==13 x₁[t]+9 x₂[t]−3 x₃[t],
  x₁[0]==−2,x₂[0]==1,x₃[0]==4},{x₁,x₂,x₃},t]]]

which yields


 {8 e^min ustmin us11 e^t + e^2 t, min us8 e^min ust + 11 e^tmin us2 e^2 t, 16 e^min ustmin us11 e^tmin use^2 t}$? >

Another way to solve this type of system of equations is to note that the equations and the initial conditions can be written in matrix form as

where

and the over dot indicates the derivative with respect to t. The solution is given by

where the square matrices [p] and [d] are solutions to [A] = [p][d][p]⁻¹. The matrices [p] and [d] are obtained from


 {p,d}=JordanDecomposition[A]

Then, the solution to the matrix form of the system of first-order equations with constant coefficients is obtained with the following program.


 A={{6,3,−2},{−4,−1,2},{13,9,−3}};
 {p,d}=JordanDecomposition[A];
 ee=p.MatrixExp[d t].Inverse[p];
 x=Flatten[ee.Transpose[{{−2,1,4}}]]//Simplify

The output of this program is identical to that obtained from DSolve, as given above.

Example 4.24

Rearrangement of a Symbolic Solution

The solution to the following fourth-order differential equation

is obtained from


 z=y[x]/.DSolve[y’’’’[x]+4 q^4 y[x]==0,y[x],x][[1]]

which displays

We shall transform this solution into one that is in terms of trigonometric and hyperbolic functions and real quantities. The interactive procedure is as follows. To obtain z in terms of trigonometric and hyperbolic functions, we use


 z1=ComplexExpand[ExpToTrig[z]]//Simplify

which gives


 C[1] Cos[(1+I) q x]+C[3] Cos[(1+I) q x]+
   (C[2]+C[4]) Cosh[(1+I) q x]+I C[1] Sin[(1+I) q x]-
   I C[3] Sin[(1+I) q x]+(−C[2]+C[4]) Sinh[(1+I) q x]

We again employ ComplexExpand to have the arguments of the trigonometric and hyperbolic function in terms of real quantities only. Then,


 z2=ComplexExpand[z1]

results in


 C[1] Cos[q x] Cosh[q x]+C[2] Cos[q x] Cosh[q x]+
  C[3] Cos[q x] Cosh[q x]+C[4] Cos[q x] Cosh[q x]-
  C[1] Cos[q x] Sinh[q x]-C[2] Cos[q x] Sinh[q x]+
  C[3] Cos[q x] Sinh[q x]+C[4] Cos[q x] Sinh[q x]+
  I (C[1] Cosh[q x] Sin[q x]-C[2] Cosh[q x] Sin[q x]-
   C[3] Cosh[q x] Sin[q x]+C[4] Cosh[q x] Sin[q x]-
   C[1] Sin[q x] Sinh[q x]+C[2] Sin[q x] Sinh[q x]-
   C[3] Sin[q x] Sinh[q x]+C[4] Sin[q x] Sinh[q x])

To simplify this result further, we collect similar terms using


 z3=Collect[z2,{Cos[q x] Cosh[q x],Cosh[q x] Sin[q x],
   Cos[q x] Sinh[q x],Sin[q x] Sinh[q x]}]

which displays


 (C[1]+C[2]+C[3]+C[4]) Cos[q x] Cosh[q x]+
   (I C[1]-I C[2]-I C[3]+I C[4]) Cosh[q x] Sin[q x]+
   (-C[1]-C[2]+C[3]+C[4]) Cos[q x] Sinh[q x]+
   (-I C[1]+I C[2]-I C[3]+I C[4]) Sin[q x] Sinh[q x]

The final step is to redefine the collection of constants. Hence


 z4=z3/.{(C[1]+C[2]+C[3]+C[4])->B₁,
  (I C[1]-I C[2]-I C[3]+I C[4])->B₂,
  (-C[1]-C[2]+C[3]+C[4])->B₃,
  (-I C[1]+I C[2]-I C[3]+I C[4])->B₄}

and the final result is


 Cos[q  x] Cosh[q  x] B₁+(Cosh[q  x] Sin[q  x] B₂+Cos[q  x] Sinh[q  x] B₃+
   Sin[q x] Sinh[q x] B₄

4.10 Solutions of Partial Differential Equations: DSolve[]

The symbolic solutions to a limited set of partial differential equations also can be obtained by using DSolve. We shall illustrate this capability with several examples. Before proceeding, it is noted that an important difference between the solutions obtained for ordinary differential equations and those obtained for partial differential equations is that the solution constants for ordinary differential equations are constants whereas those for partial differential equations may be functions of one or more of the independent variables.

Example 4.25

Solution to Partial Differential Equation #1

Consider the partial differential equation

Its solution is obtained with the following statement


 sol=u[x,y]/.DSolve[x D[u[x,y],x]-y D[u[x,y],y]+
   y^2 u[x,y]==y^2,u[x,y],{x,y}][[1]]

which gives

It is seen that the C[1] is a function of the product of x and y.

To obtain the value of sol when x = 2 and y = 3, we use


 sol/.{x->2.,y->3.}

to obtain


 1+90.0171 C[1][6.]

Example 4.26

Solution to Partial Differential Equation #2

Consider the partial differential equation

with the boundary condition v(x,0) = sinx. Its solution is obtained with


 sol2=v[x,t]/.Simplify[DSolve[{D[v[x,t],x]+D[v[x,t],t]+
   σ v[x,t]==0,v[x,0]==Sin[x]},v[x,t],{x,t}]][[1]]

which yields


min use^min usσtSin[tmin usx]

Example 4.27

Solution to Partial Differential Equation #3

Consider the partial differential equation

Its solution is obtained with


 sol3= w[x,y]/.DSolve[D[w[x,y],x,x]+5 D[w[x,y],x,y]+
  6 D[w[x,y],y,y]==0,w[x,y],{x,y}][[1]]

which gives


 C[1][−2 x+y]+C[2][−3 x+y]

4.11 Laplace Transform: LaplaceTransform[] and InverseLaplaceTransform[]

The Laplace transform is obtained from


 LaplaceTransform[expr,t,s,Assumptions->{…}]

where s is the Laplace transform parameter, t is the independent variable that is being transformed, expr is an expression that is a function of t or is a constant, and Assumptions is used to place restrictions on any parameters appearing in expr. The initial conditions are dealt with by using Simplify and its Assumptions option.

Before proceeding, it is again noted that the definition of the unit step function u(x) as given by UnitStep[x] is used in Mathematica to represent a piecewise function such that u(x) = 0 for x < 0 and u(x) = 1 for x ≥ 0. Note that this definition includes 0. This is different from the Heaviside Theta function, which is given by HeavisideTheta. As mentioned in Example 4.20, this function is defined as 0 when its argument is < 0 and equal to 1 when its argument is > 0. It is not defined for an argument equal to 0.

We now consider the following example.

Example 4.28

Laplace Transform of a Half Sine Wave

We shall determine the Laplace transform of the half sine wave

which is obtained from


 lt=Simplify[LaplaceTransform[(1-UnitStep[t-π/a]) Sin[a t],
  t,s,Assumptions->a>0]]

Its execution gives

The inverse Laplace transform is obtained from


 ilt=InverseLaplaceTransform[lt,s,t]

which displays

To verify that ilt is correct, we use Simplify with two different assumptions as follows


 Simplify[ilt,Assumptions->t<π/a]
 Simplify[ilt,Assumptions->t>π/a]

which gives, respectively,


 Sin[a t]
 0

Example 4.29

Laplace Transform Solution of an Inhomogeneous Differential Equation #1

We continue with Example 4.28 in the following manner. Consider the ordinary differential equation

(a)

with y(0) = y′(0) = 0. Because of the way the general solution from LaplaceTransform is presented, we shall obtain the Laplace transform of this equation and its inverse one expression at a time. The Laplace transform of Eq. (a) is obtained from


 ltde=Simplify[LaplaceTransform[y’’[t]+y[t]==
   (1-UnitStep[t-π/a]) Sin[a t],t,s,Assumptions->a>0],
   Assumptions->{a>0,y[0]==0,y’[0]==0}]

which displays

The unknown quantity is LaplaceTransform[y[t],t,s,Assumptions->a>0]. This is a little awkward to use in this form. Therefore, to improve readability, we shall use the transformation rule and redefine it (arbitrarily) as laps. Then, we start over and use the modified instruction


 ltde=Simplify[LaplaceTransform[y’’[t]+y[t]==
   (1-UnitStep[t-π/a]) Sin[a t],t,s,Assumptions->a>0]/.
   LaplaceTransform[y[t],t,s,Assumptions->a>0]->laps,
   Assumptions->{a>0,y[0]==0,y’[0]==0}]

which displays

Next, we solve for the Laplace transform of y(t), which is now represented by the quantity laps, by using Solve as follows


 sols=laps/.Simplify[Solve[ltde,laps]][[1]]

which results in

The inverse Laplace transform is determined from


 ytt=InverseLaplaceTransform[sols,s,t]

which yields

It is noted that when a → 1


 Limit[ytt,a->1]

we obtain

To plot ytt with, say, a = 2, we use


 a=2.;
 Plot[ytt,{t,0,30},AxesLabel->{"t","y(t)"}]

which results in Figure 4.4.

Example 4.30

Response of a Two Degrees-of-Freedom System

The nondimensional governing equations of motion of a two degrees-of-freedom system are [1]

where x_j is the displacement of mass m_j as a function of the nondimensional time τ. The quantities ω_r, m_r, and ζ_j are nondimensional parameters, k₁ is the spring stiffness, and f_j is the force applied to mass m_j.

We shall solve this system of ordinary differential equations using the Laplace transform. It is assumed that the initial conditions are zero; that is, x₁(0) = dx₁(0)/dτ = x₂(0) = dx₂(0)/dτ = 0, and that f₁(t) = 0. In addition, to improve readability of the results we shall replace the system’s output variable indicating the Laplace transform of the function by a shorter designation using the transformation rule. In particular, we shall use the following transformations


 LaplaceTransform[x1[t],t,s]->x1s
 LaplaceTransform[x2[t],t,s]->x2s
 LaplaceTransform[f1[t],t,s]->f1s
 LaplaceTransform[f2[t],t,s]->f2s

Then the Laplace transform of the governing equations are obtained with


 lt1 = Simplify[LaplaceTransform[
 x₁’’[t]+(2 ζ₁+2 ζ₂ m_r ω_r) x₁’[t]+(1+m_r ω_r^2) x₁ [t]
 −2 ζ₂ m_r ω_r x₂’[t]-m_r ω_r^2 x₂[t]==f₁[t]/k₁,t,s]/.
 {LaplaceTransform[x₁[t],t,s]->sx₁,
 LaplaceTransform[x₂[t],t,s]->sx₂,
 LaplaceTransform[f₁[t],t,s]->sf₁},
 Assumptions->{x₁’[0]==0,x₁[0]==0,x₂[0]==0,sf₁==0}]
lt2=Simplify[LaplaceTransform[
 x₂’’[t]+2 ζ₂ ω_r x₂’[t]+ω_r^2 x₂[t]−2 ζ₂ ω_r x₁’[t]
 -ω_r^2 x₁[t]==f₂[t]/(k₁ m_r),t,s]/.
 {LaplaceTransform[x₁[t],t,s]->sx₁,
 LaplaceTransform[x₂[t],t,s]->sx₂,
 LaplaceTransform[f₂[t],t,s]->sf₂},
 Assumptions->{x₂’[0]==0,x₂[0]==0,x₁[0]==0}]

which gives

We now solve for sx₁ and sx₂ using


 sols2=Flatten[{sx₁,sx₂}/.Simplify[Solve[lt1&&lt2,
{x1s,x2s}]]]

which results in the two-element list

The symbolic inverse Laplace transform of these results is not very useful. Therefore, we shall assign numerical values to the constants and assume that the force on m₂ is an impulse force; that is, f₂(t) = δ(t). The parameters are: _r = 0.85, m_r = 0.45, k₁ = 1, and ζ₁ = ζ₂ = 0.15. Then the plotted inverse Laplace transforms are obtained with


 f2s=LaplaceTransform[DiracDelta[t],t,s];
 m_r=0.45; ω_r=0.85; ζ₁=0.15; ζ₂=0.15; k₁=1;
 ilt1=InverseLaplaceTransform[sols2[[1]],s,t];
 ilt2=InverseLaplaceTransform[sols2[[2]],s,t];
 Plot[{Chop[ilt1],Chop[ilt2]},{t,0,50},PlotRange->All,
  AxesLabel->{"τ","x₁(τ), x₂(τ)"}]

which results in Figure 4.5.

Example 4.31

Laplace Transform Solution of an Inhomogeneous Differential Equation #2

Consider the following differential equation

where

We shall solve this ordinary differential equation using the Laplace transform. It is assumed that the initial conditions are zero; that is, r(0) = dr(0)/dt = 0. In addition, to improve readability of the results we shall replace the system’s output variable indicating the Laplace transform of the function by a shorter designation using the transformation rule. Then the statements to obtain the Laplace transform of the governing equations are


 rs=Simplify[LaplaceTransform[r’’[t]+2 p r’[t]+r[t]==
  Piecewise[{{1,0<=t<t1},{−1,t1<=t<2 t1}}],t,s]/.
  LaplaceTransform[r[t],t,s]−>lts,
  Assumptions−>t1>0&&r[0]==0&&r’[0]==0];
  rts=lts/.Solve[rs,lts][[1]]

The output of this program is

To determine the solution in the time domain, we take the inverse Laplace transform using


 rt=Simplify[InverseLaplaceTransform[rts,s,t]];

The result is rather lengthy and has been omitted. However, it can be plotted. Therefore, we shall assume that t₁ = 0.5 and p = 0.15 and the plotted results, which are shown in Figure 4.6, are obtained with


 t1=0.5; p=0.15;
 Plot[Chop[rt],{t,0,25},AxesLabel−>{"t","r(t)"}]

4.12 Functions Introduced in Chapter 4

Listed in Table 4.5 are the functions introduced in this chapter that do not appear in Table 4.1.

Table 4.5 Summary of additional commands introduced in Chapter 4

Command	Usage
Coefficient	Gives the coefficient of a polynomial or power series
CoefficientList	Gives a list of the coefficients of a polynomial or a power series
D	Determines the ordinary or partial derivative or multiple derivative of a symbolic function
DiracDelta	Represents the delta function
Dsolve	Attempts to obtain a symbolic solution to an ordinary or partial differential equation
EulerEquations	Obtains the Euler–Lagrange equations [Requires loading of the Variational Methods Package]
HeavisideTheta	Represents the Heaviside theta function of argument x: when x < 0, output is 0; when x > 0, output is 1; undefined at x = 0
Integrate	Attempts to obtain a symbolic solution of an indefinite or definite integral
InverseLaplaceTransform	Attempts to determine the inverse Laplace transform
JordanDecomposition	Obtains [p] and [d] as the solution to [A] = [p][d][p]⁻¹ for square matrices [A], [p], and [d] when [A] is given
LaplaceTransform	Attempts to determine the Laplace transform of an expression or an ordinary differential equation
Limit	Finds the limiting value of an expression
Normal	Converts some special forms of results to a normal form
Piecewise	Represents a piecewise function over the regions specified in its argument
Residue	Finds the residue of an expression at one of its singularities
Series	Gives the power series expansion of an expression or of a common function
Solve	Attempts to find a symbolic solution to an equation or system of equations
UnitStep	Represents the unit step function of argument x: when x < 0, output is 0; when x ≥ 0, output is 1

References

B. Balachandran and E. B. Magrab, Vibrations, 2nd edn, Cengage Learning, Toronto, Ontario, 2009, p. 472.
S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd edn, McGraw-Hill, New York, 1959, pp. 298–300.

Exercises

Table 4.1

4.1 The command TrigExpand and TrigReduce also work with hyperbolic functions. Consider the expression

Show that this can be reduced to
4.2 Find an expression for the real part of w in terms of trigonometric and hyperbolic functions when

Section 4.3

4.3 Solve for x when
4.4 Given the function

The value of x = x_max is that value at which the maximum value of f(x) occurs. It can be determined from the solution to df(x)/dx = 0 provided that d²f(x_max)/dx² < 0. Determine x_max.
4.5 Given the following system of equations

Determine expressions for i_k.
4.6 Given the following set of equations

Solve for x_j and then use the results in the following equation to obtain a cubic equation in h from which h can be determined

Do not solve for h.
4.7 If

then determine an expression for α, α > 0, that satisfies f_m(α) = f_m₊₁(α)
4.8 The nondimensional buckling load α of bar under a distributed axial load and resting on an elastic foundation can be determined from the determinant of the coefficients a_j of

where γ is a nondimensional quantity that is proportional to the spring modulus of the elastic foundation.
1. Find the symbolic solution for α when K = 1 and when K = 3.
2. Find the numerical value of α when K = 1, 3, and 5 and γ = 0 and 0.4. Present the results in a table with the appropriate labels.
4.9 Given the following determinant

where Δ = Ω² is a nondimensional frequency coefficient. Determine the smallest positive value of Ω when k_a = π/4, ν = 0.3, and n = 0, 1, …, 10. Display the results in tabular form.

Section 4.4

4.10 Find the following limits
4.11 Find the limit of the following expression when n is an integer

Section 4.5

4.12 Obtain the coefficients of a five-term expansion of the following expression about = 0.
4.13 Obtain the coefficients of the polynomial in γ that result from expanding the determinant
4.14 Find the sum of the squares of the coefficients of the first five terms of the series expansion around x = 0

Section 4.7

4.15 Using the expression for the curvature of a parametric curve given in Example 4.10, show that when

the curvature is given by

Section 4.8

4.16 Evaluate the following integral when b > 0 and b is real.
4.17 Evaluate the following integral, which is proportional to the total radiation by a black body at all frequencies.

Section 4.9

4.18 The nondimensional transverse displacement y of a uniformly loaded annular circular plate with Poisson’s ratio ν is given by

Determine an expression for y when the boundary conditions at the inner boundary ξ = α, 0 < α < 1, are

and those at the outer boundary ξ = 1 are

Plot the results for y(ξ), y′(ξ), y′′(ξ), and y′′′(ξ) when α = 0.2 and ν = 0.3 using the grid plot commands of Example 4.20.
4.19 Determine the general solution to the following equation

where w = w(r), q_o is a constant, and m = 0, 1, …, 5.
4.20 Consider a solid circular plate whose thickness varies as

where β is a constant. The governing equation for this type of plate when it is subjected to a nondimensional uniform load of magnitude p over the entire plate is given by [2, p. 298–300]

and w is the transverse displacement.
1. Verify that a particular solution to this equation is
2. Determine the solution to the homogeneous equation; that is, when p = 0.
3. If the solution is to remain finite at η = 0, then take the solution obtained in (b) and determine which part of the solution remains finite at η = 0.
4.21 In terms of nondimensional quantities, the governing equation of an annular circular plate whose thickness varies linearly and is subjected to a uniform load Q_o distributed over the entire plate and a line load along the inner boundary of magnitude P_o is given by [2, p. 303–305]

where ν is the Poisson’s ratio. Determine the solution to this equation when ν ≠ 1/3 and when ν = 1/3.
4.22 Given the following differential equation

where 1. Assume that an approximate solution to this equation is
1. Substitute this solution into the differential equation and generate a list of three equations of the form
2. Obtain a solution to these three, coupled equations.

Section 4.11

4.23 Find the inverse Laplace transform of

and plot the result.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.

Table of Contents for An Engineer’s Guide to Mathematica®

Create new playlist

Sign In

Sign Up

4.1 Introduction

4.2 Assumption Options

4.3 Solutions of Equations: Solve[]

4.4 Limits: Limit[]

4.5 Power Series: Series[], Coefficient[], and CoefficientList[]

4.6 Optimization: Maximize[]/Minimize[]

4.7 Differentiation: D[]

4.8 Integration: Integrate[]

4.9 Solutions of Ordinary Differential Equations: DSolve[]

4.10 Solutions of Partial Differential Equations: DSolve[]

4.11 Laplace Transform: LaplaceTransform[] and InverseLaplaceTransform[]

4.12 Functions Introduced in Chapter 4

References

Exercises

Table of Contents for
An Engineer’s Guide to Mathematica®