Instructions for Math 251, Applied Differential Equations. Only selected

problems from the list below are collected for grading. There is a separate schedule, updated weekly, for these collections. All work is collected in class and returned in class, with the exception that work can be submitted under the door at JWB 113 any time on the due date.

Late Work. The only work submitted in class is work that is on time. Catch-ups are slipped under the door at JWB 113. There is a two–day grace period for receipt of work. After two days it is recorded as a zero. The lowest 5 scores are dropped, to eliminate the need to catch up when disaster strikes you. Once you get 5 zero scores, then it is time to call 581-6879. I’ll schedule a summit conference to solve the problem.

Working Together. A group may consult with other groups to check on answers or methods. A group is to submit one solution, not multiple solutions. Group size is 1 to 3.

Handwritten Solutions. All problems are to be solved at a desk with pencil and paper. The rules outlined in the syllabus apply to these presentations. In particular, they are not scratch paper solutions or just answers. The solutions represent considerable thought and effort, in short, your very best work.

Maple may be used as an assist when appropriate to check steps or to do details that are subject to human error (Less than 5% of the course has anything to do with maple). In particular, maple is not to blame for wrong answers, you are to blame. Ditto for hand calculators and other assists to computation or graphing.

Partial Credit. Wrong answers receive partial credit based upon what is written on the paper, the steps shown and the logic exhibited. You will be taught by example what is expected in a solution. Generally, the textbook is too brief, and it is a bad example of how to write up a solution. References by theorem name and page are important. Sentences that explain mathematical steps are essential to written solutions.

Format. Submit solutions in one-sided format, stapled in the upper left corner. A problem starts on a new page, with the statment first. Observe the usual rules for solution: state the problem, solve it, and give all details, especially page references to the textbook. The Bradley Gigi format described in the syllabus and engineering paper are preferred. But there is no strict rule about this.

Graphs. A graphic from your hand calculator is supposed to be a hand–drawn replica. Maple graphs are expected to be presented in compressed format, about 2×2 inches, preferably pasted-up onto your engineering paper, or replicated by hand. A screen print may be used to get the size right (see the problem notes) or else trace the graph from the screen or draw a likeness by hand.

Maple and Calculator Results. Generally, the solution to a problem is completely devoid of the mention of maple, even though it may have been used at some point in the solution, just to check answers. This is normal. Don’t force maple into the solution. Maple session output is for you, not for the report. Keep the report handwritten and write a note that you checked it with maple or derived the result with maple commands. Pencil in the essence of the maple idea, not a log script. The same remarks apply to graphing calculators.

Math 251 Take–Home Exams

Submit solutions only for the problems listed in the schedule. The schedule is updated weekly.

Separable equations. Obtain the implicit form of the solution, and if possible, also the explicit form of the solution.

Problem 1. x′ = −1 + exp(−2x), x(0) = 1. Problem 2. x′ = t/(x^{2} + 1), x(0) = 1.

Problem 3. y′ = \frac{xy + 3x − y − 3}{xy − 2x + 2y − 4}, y(0) = 1.

Problem 4. dP/dt = P (1 − ln(P)), P (0) = 1. Problem 5. dT/dt = −0.7(T − 75), T (0) = 180. Problem 6. (1 − t/10)dI/dt = 4 − 0.4I, I(0) = 0. Problem 7. y′ = y^{2} + 5y + 4, y(0) = 10. Problem 8. yy′ sec^{2}(x) = y^{2} + 1, y(0) = 1.

Linear equations. Solve for the solution subject to the given initial condition. Undefined variables are constants.

Problem 9. (1 + x^{2})y′ + y = 0, y(0) = 1 Problem 10. x′ − x = t, x(0) = 1 Problem 11. x′ + x = sin(t), x(0) = 0 Problem 12. x′ + (1/t)x = t, x(1) = 1 Problem 13. dU/dt = A − BU, U(0) = 0 Problem 14. m(dv/dt) = mg − kv, v(0) = 0 Problem 15. dy/dx + (1/x)y = 3, y(1) = 0 Problem 16. x(dw/dx) − 4w = x^{6} exp(x), w(1) = 0

Exact differential equations. Test for exactness. If exact, then solve for the implicit form of the solution. If possible, obtain also the explicit form of the solution.

Problem 17. (5x + 4y)dx + (4x − 8y^{3})dy = 0 Problem 18. (tan x − sin x sin y) + cos x cos y \frac{dy}{dx} = 0 Problem 19. (x ln x − exp(−xy))\frac{dy}{dx} + (x^{−1} + y ln x) = 0

Problem 20. (1 − \frac{3}{y} + x)\frac{dy}{dx} + (1 − \frac{3}{x} + y) = 0

Classification of first order differential equations. Classify as separable, linear, exact

and solve:

Problem 21. x′ = −1 + exp(−x), x(0) = 1 Problem 22. x′ = (t + 1)/(x^{4} + 1), x(0) = 1 Problem 23. x′ = x + t, x(0) = 1 Problem 24. x′ = −x + 4 sin(t), x(0) = 0 Problem 25. x′ = (−4/t)x + t^{4}, x(1) = 0 Problem 26. 2xy + (1 + x^{2})y′ = 0, y(0) = 1 Problem 27. e^{y} + (xe^{y} − 1)y′ = 0, y(5) = 0 Problem 28. y′ = 1 − y^{2}, y(0) = 0

Applications of first order equations.

Problem 29. (LR Circuit)Solve the equation LI′(t) + RI(t) = E(t), I(0) = 0, given R = 10,

L = 2, E(t) = 200. Problem 30. (RC Circuit) Solve the equation RQ′(t) + (1/C)Q(t) = E(t), Q(0) = 0, given

R = 1, C = 0.01, E(t) = sin(60t). Problem 31. (Phosphate Depletion)The amount P (t) in grams of phosphate in a certain

human body satisfies the equation \frac{dP}{dt} = −\frac{1}{3} + \frac{1}{12} cos(π(t − 6)/24), P (0) = 600. Find the amount of phosphate P (24) in the body after 24 hours.

Problem 32. (Radioactive Waste)Strontium 90 has a half-life of 25 years. Find the percentage of radioactive material left in a sample after a decay period of 32 years.

Fundamental Theorem of Calculus. Solve the differential equation exactly and also numerically on [0, 1]. Present a table for the numerical solution and graph the solution from the table values.

Problem 33. (Trapezoidal Rule)y′ = x exp(−x), y(0) = 0, Trapezoidal Rule step 0.05. Problem 34. (Simpson Rule)y′ = x ln(1 + x), y(0) = 1, Simpson Rule step 0.05. Problem 35. (Trapezoidal Rule)y′ = x^{5}(1 − x^{3}), y(0) = 0, Trapezoidal Rule step 0.05.

Problem 36. (Simpson Rule)y′ = π sin^{2}(πx) cos(πx), y(0) = 0, Simpson Rule step 0.05.

Applications.

Problem 37. (Cable Throw)A cable 1000 feet long weighs 1.5 pounds per linear foot. On one end is attached a metal ball weighing 30 pounds. The ball is thrown vertically upward with velocity v feet per second, but the weight of the cable affects the upward progress x(t) of the ball. It is known that

\frac{d}{dt} ([30 + 1.5x(t)]x′(t)) = −32(30 + 1.5x(t)),

x(0) = 0, x′(0) = v.

Let y(v) be the maximum height attained by the ball for initial velocity v.

Without solving the differential equation show that y and v satisfy the equation

32(30 + 1.5y)^{3} = (45v)^{2} + 32(30)^{3}

Determine v = v_{0} for which y = 900 and plot y versus v on 0 ≤ v ≤ v_{0}.

Problem 38. (Rectilinear Motion with Drag Force)A rectilinear motion problem is modeled with drag force proportional to the velocity squared, giving the specific model for the velocity v(t):

\{ \begin{matrix}v′ + (24/32)v^{2} = 32, \\ v(0) = 0. \end{matrix}

Solve exactly for v(t) in the problem v′ = 32 − (24/32)v^{2}, v(0) = 0, 0 ≤ t ≤ 2.0 and Graph v(t) versus t on 0 ≤ t ≤ 2.0.

Constant Coefficient Equations. Apply Euler’s method to determine the Euler auxiliary equation and then find the general solution y = c_{1}y_{1} + c_{2}y_{2} from the recipe for constant equations.

Problem 39. y′′ = 0

Problem 40. y′′ + y′ = 0

Problem 41. y′′ + y = 0

Problem 42. y′′ + 3y′ + 2y = 0

Problem 43. y′′ + 2y′ + y = 0

Problem 44. y′′ + y′ + y = 0

Problem 45. y′′ − 9y = 0

Problem 46. 9y′′ + y = 0 Euler Differential Equations. Use the recipe for Euler-Cauchy equations to calculate a solution basis y_{1}, y_{2} of the given Euler differential equation.

Problem 47. x^{2}y′′ = 0. Problem 48. x^{2}y′′ + xy′ = 0. Problem 49. x^{2}y′′ − xy′ = 0. Problem 50. x^{2}y′′ − 2xy′ + 2y = 0. Problem 51. 4x^{2}y′′ + y = 0. Problem 52. x^{2}y′′ − xy′ + 2y = 0. Problem 53. 4x^{2}y′′ − 4xy′ + 3y = 0. Problem 54. x^{2}y′′ − 5xy′ + 9y = 0.

Problem 55. (Wronskian) Given y_{1} = x^{3}, y_{2} = x^{4} and the Euler differential equation x^{2}y′′ −

6xy′ + 12y = 0, verify that y_{1}, y_{2} are linearly independent solutions on (0, ∞). Problem 56. (Wronskian Test) Given y_{1} = x^{3}, y_{2} = |x^{3}| and the Euler differential equation

x^{2}y′′ − 4xy′ + 6y = 0, verify that y_{1}, y_{2} are linearly independent solutions on (−∞, ∞) with identically zero Wronskian. Does this violate the Wronskian Test?

Problem 57. (Superposition) Let x_{1}(t) be the response to external force g_{1}(t) in the pendulum equation x′′(t) + a^{2}x(t) = g_{1}(t), x(0) = x′(0) = 0. Suppose x_{2}(t) is the response for force g_{2}(t), x′′(t) + a^{2}x(t) = g_{2}(t), x(0) = c, x′(0) = d. Verify that the sum of the responses

y(t) = x_{1}(t) + x_{2}(t) satisfies y′′ + a^{2}y = g_{1}(t) + g_{2}(t), y(0) = c, y′(0) = d.

Problem 58. (Reduction of Order) Given xy′′ + y′ = 0 and a solution y_{1}(x) = ln(x), find a basis

of solutions on (0, ∞) by the method of reduction of order. Problem 59. (Reduction of Order) Given (1 − x^{2})y′′ − xy′ = 0 and the solution y_{1}(x) = 1,

determine a basis of solutions by the method of reduction of order.

Forced Equations. For each equation determine the characteristic equation and the general solution c_{1}y_{1} + c_{2}y_{2} of the unforced equation y′′ + ay′ + by = 0 and a particular solution y_{p} of the forced equation y′′ + ay′ + by = g(x). Apply either the method of undetermined coefficients or the method of variation of parameters as indicated.

Problem 60. y′′ − 2y′ − 3y = x − 4 + 2xe^{2x} by undetermined coefficients. Problem 61. y′′ + 3y′ + 2y = e^{x} by both methods. Problem 62. y′′ − 8y′ + 25y = x by the easiest method. Problem 63. y′′ − 4y = e^{2x} by variation of parameters.

Problem 64. y′′ + y = 8 sin^{2}(x) by either method.

Problem 65. y′′ + y = sin(x) cos(x) by either method.

Second Order Applications.

Problem 66. Solve the Cantilever Equation EIy′′ = \frac{w(L − x)^{2}}{2} , y(0) = y′(0) = 0.

Problem 67. Solve the Overdamped Spring Equation x′′(t) + 3x′(t) + 2x(t) = 0.

Problem 68. Solve the Critically Damped Spring Equation x′′(t) + 16x′(t) + 64x(t) = 0. Explain in words why the classification is changed if the constant 16 is replaced by 16.1.

Problem 69. Solve the Underdamped Spring Equation x′′(t) + 2x′(t) + 2x(t) = 0.

Problem 70. Solve the forced, damped LRC circuit equation (0.2)x′′ + (1.2)x′ + 2x = 4 sin(4t), x(0) = 0, x′(0) = 0 and identify the steady state and transient solutions.

Problem 71. Solve the forced, damped torsional motion equation (for a twisting shaft) Iθ′′ + cθ′ + kθ = T (t), θ(0) = 0, θ′(0) = 0 for I = 2, c = 2, k = 3, T (t) = sin(t).

Problem 72. (RLC Circuit)The RLC circuit equation with initial conditions at t = 0 is (′ = d/dt):

LQ′′ + RQ′ + \frac{1}{C} Q = E(t), Q(0) = q_{0}, Q′(0) = q_{1}.

Assume the specific form

E(t) = e_{1} cos(wt) + e_{2} sin(wt)

Solve the differential equation exactly and plot the solution on [0, T] for the following data sets:

\begin{matrix}Set & T & L & R & C & e_{1} & e_{2} & w & q_{0} & q_{1} \\ 1 & 5 & 1 & 3 & 1 & 1 & 1 & 10 & 0 & 1 \\ 2 & 5 & 1 & 2 & 1 & 1 & -1 & 10 & 0 & 1 \\ 3 & 5 & 1 & 1 & 1 & -1 & 1 & 10 & 1 & 0 \end{matrix}

Mechanical Resonance. This study of resonance considers the forced harmonic oscillator with zero initial states:

\{ \begin{matrix}x′′ + a^{2}x = g(t), \\ x(0) = x′(0) = 0, \end{matrix}

Problem 73. (Convolution formula)Verify that x is given by the variation of parameters

solution

x(t) = ∫_{0}^{t} \frac{1}{a} sin(a(t − r))g(r)dr.

The integral on the right has the form ∫_{0}^{t} f(t − r)g(r)dr, which is called the convolution of the functions f and g.

Problem 74. (Resonance) According to the theory of mechanical resonance the forcing term g(t) = 2a^{2} sin(at) will cause the amplitudes of x(t) to grow without bound as t → ∞. Illustrate the theory by verifying that x(t) = sin(at) − at cos(at) is the solution (it has amplitude \sqrt{1 + a^{2}t^{2}}).

Problem 75. (Wine glass experiment)Suppose x(t) approximates the deflection from equilibrium of the rim of a wine glass subjected to sound waves from a loudspeaker. The Hooke’s law term a^{2}x(t) would be derived from elastic properties of the wine glass. The forcing term g(t) would be derived from the sound waves emitted by the loudspeaker. Explain in a short paragraph how the theory of resonance applies to design an experiment to break the wine glass.

Higher Order Equation Solution Basis. Find a solution basis for the given differential equation. First try to factor the characteristic polynomial into linear and quadratic factors. If this is possible then an exact basis should be given. If not then compute and use approximate roots.

Problem 76. x′′′ + 3x′′ + 4x′ + 2x = 0.

Problem 77. 100x′′′ − 334x′′ + 668x′ − 536x = 0.

Problem 78. 100x′′′ − 334x′′ + 668x′ − 535x = 0.

Problem 79. 9x^{(4)} − 42x^{(3)} + 22x^{(2)} + 126x^{(1)} − 147x = 0.

Initial Value Problem. Calculate the solution of the initial value problem with the prescribed initial values.

Problem 80. x′′′ − 4x′′ − 80x′ − 192x = 0, x(0) = 1, x′(0) = 0, x′′(0) = 0.

Problem 81. 4x′′′ − 3x′′ − 18x′ + 17x = 0, x(0) = 0, x′(0) = 1, x′′(0) = 0.

Problem 82. x′′′ + 6x′′ + 11x′ + 6x = 0, x(0) = 1, x′(0) = 1, x′′(0) = 1.

Problem 83. x^{iv} + 2x′′′ + 2x′′ + 2x′ + x = 0. x(0) = 1, x′(0) = 0, x′′(0) = 0, x′′′(0) = 1.

Power Series Method. Substitute the trial power series solution

x(t) = ∑_{n=0}^{∞} c_{n}t^{n}

into the given initial value problem. Find a general recursion formula for c_{n} in terms of c_{0} and c_{1}, then solve for the coefficients c_{0}, c_{1}, ⋯, c_{7} using the given initial conditions.

Problem 84. x′′ + tx = 0, x(0) = 0, x′(0) = 1.

Problem 85. x′′ + 2tx′ − x = 0, x(0) = 0, x′(0) = 1.

Problem 86. x′′ + e^{t}x = 0, x(0) = 1, x′(0) = 0. Problem 87. x′′ + t^{2}x = 0, x(0) = 0, x′(0) = 1.

Legendre Polynomials. The Legendre equation

(1 − x^{2})y′′ − 2xy′ + μy = 0, −1 ≤ x ≤ 1,

is known to have a non-trivial solution that is finite when x = ±1 if and only if μ = n(n + 1) where n = 0, 1, 2, · · ·. That solution y(x) = P_{n}(x) is a polynomial of exact degree n and is odd (resp. even) if n is odd (resp. even). It is customary to normalize P_{n}(x) by the condition

P_{n}(1) = 1, n = 0, 1, 2, 2, · · · .

With this convention the polynomials are uniquely defined. Power series methods can be used to show that they are related by the Legendre polynomial recurrence relation

(n + 1)P_{n+1}(x) = (2n + 1)xP_{n}(x) − nP_{n−1}(x), n ≥ 1.

Problem 88. Use power series methods to directly verify that P_{0}(x) = 1 and P_{1}(x) = x. Problem 89. Use the Legendre polynomial recurrence relation to calculate P_{2}(x) through P_{8}(x). Problem 90. Use the Legendre polynomial recurrence relation and mathematical induction to

verify that P_{n}(1) = 1 for all n = 0, 1, 2, · · ·. Problem 91. Produce graphs of P_{0}(x), P_{1}(x), …, P_{8}(x) on −1 ≤ x ≤ 1.

Computation of Transforms. Compute using the theory of the transform and transform tables.

Problem 92. 3t + 2 Problem 93. t^{2} − 4t + 1.1 Problem 94. (t − 1)(t + 1)(t − 2) Problem 95. 5e^{1−t} Problem 96. e^{−t/3} + sinh(−t/3) Problem 97. (t^{2} − t)e^{−t} Problem 98. e^{t} sin(t) − e^{2t} cos(4t) Problem 99. t^{2} cos(3t) Problem 100. te^{−t} sin(2t) Problem 101. 2te^{−3t} cosh(t)

Computation of Transforms. Compute using the theory of the transform, the theory for Heaviside’s unit step function H and transform tables.

Problem 102. sin^{2}(3t) Problem 103. sin(t)/t Problem 104. (1 − cosh(2t))/t Problem 105. e^{t} + cos(t − π)H(t − π) Problem 106. H(t) + H(t − 1) + H(t − 2) + H(t − 3) Problem 107. | sin(2t)| Problem 108. The 2-periodic extension of H(t) + H(t − 1) on [0, 2] to the whole real line. Problem 109. The 2-periodic extension of tH(t) + (1 − 2t)H(t − 1) on [0, 2] to the whole real

line.

Computation of Inverse Transforms. Write each function of s as the Laplace transform L(f(t)) for some function f(t). Cite the use of basic theorems.

Problem 110. \frac{3}{s^{2}}
Problem 111. \frac{s + 1}{s^{2} + 1}
Problem 112. \frac{1}{(s − 2)^{3}}
Problem 113. \frac{3s + 2}{s^{2} + 6s + 25}
Problem 114. \frac{s − 1}{(s − 1)^{2} − 4}
Problem 115. \frac{s + 1}{(s − 1)(s − 2)}
Problem 116. ln (1 + \frac{4}{s^{2}} )

Problem 117. arctan(1/s) Problem 118. \frac{1}{s} arctan(1/s)

Problem 119. \frac{exp(−2s)}{s}

Computation of Inverse Transforms. Write each function of s as the Laplace transform L(f(t)) for some function f(t). Cite the use of basic theorems.

Problem 120. \frac{2e^{−s}}{s^{2} + 16}

Problem 121. \frac{1}{s}(1 − exp(−2s))^{2}

Problem 122. \frac{1}{s} tanh(s)

Problem 123. \frac{exp(−2s)}{s^{2} − 1}
Problem 124. \frac{1 − exp(−4s)}{s^{2}}
Problem 125. \frac{s(1 − exp(−2πs))}{s^{2} + 1}

Problem 126. \frac{1}{(s^{2} + 1)^{3}} + s^{−5} exp(−10s)

Problem 127. \frac{1}{s^{4}(s^{2} + 1)}

Solution of Differential Equations. Compute the transformed equation for L(y) and also y(t) by Laplace transform methods.

Problem 128. y′′ + 4y = 0, y(0) = 0, y′(0) = 2 Problem 129. y′′ − 4y = 0, y(0) = 1, y′(0) = 1 Problem 130. y′′ + 3y′ + 2y = 0, y(0) = 0, y′(0) = 1 Problem 131. y′′ + 3y′ + 2y = exp(−t), y(0) = y′(0) = 0 Problem 132. y′′ − 5y′ + 4y = exp(3t), y(0) = 1, y′(0) = −1 Problem 133. y′′ + 2y′ + 2y = t + sin(t), y(0) = y′(0) = 0

Differential Equation Methods. Find the transform equation and compute y(t) by Laplace transform methods.

Problem 134. L(y) = Y (s), Y′ + (1/(s − 1))Y = 0, Y (0) = 2 Problem 135. ty′′ − ty′ − y = t, y(0) = 0, y′(0) = 1 Problem 136. y(t) = Q(t), E(t) = 10 sin(10t), LQ′′ + RQ′ + (1/C)Q = E(t), Q(0) = Q′(0) = 0

L = 1, R = 6, C = 1/9 Problem 137. y′′′ + y = 0, y(0) = y′(0) = 0, y′′(0) = 1 Problem 138. y^{(4)} − y = 0, y(0) = 0, y′(0) = y′′(0) = 0, y′′′(0) = −1 Problem 139. y^{(4)} − y = cos(t), y(0) = y′′(0) = 1, y′(0) = y′′′(0) = 0

Solution of Special Equations. Find the transform L(y) and compute y(t) using Laplace transform methods for the Heaviside function and the convolution integral.

Problem 140. \frac{d}{ds}L(y) = \frac{∫_{0}^{π} e^{−st} sin(t)dt}{1 − exp(−πs)}

Problem 141. \frac{d^{2}}{ds^{2}} L(y) = \frac{2 e^{−2 s−2}}{s + 1} + \frac{e^{−2 s−2}}{(s + 1)^{2}}

Problem 142. y′′ + y = H(t − 10),y(0) = y′(0) = 0 Problem 143. y′′ + 10y′ + 50y = E(t), y(0) = y′(0) = 0, E(t) = 25t(1 − H(t − 4)) + 50H(t − 4) Problem 144. y(t) = exp(−t) − ∫_{0}^{t} sin(t − r)y(r)dr

Problem 145. y(t) = −t + ∫_{0}^{t} \frac{(t − r)^{3}}{6} y(r)dr

Partial Fraction Applications. Solve for y(t) in the following equations using the partial fraction methods of Oliver Heaviside and properties of Laplace transforms.

Problem 146. L(Q)[s^{2} + 6s + 9] = \frac{s + 1}{s^{2} + 1}

Problem 147. L(y) = \frac{s}{(s^{2} + 1)(s + 1)^{2}}

Problem 148. L(y) = \frac{2}{(s + 1)(s − 1)(s + 2)}

Problem 149. L(y) = \frac{exp(−4s)}{(s − 1)(s + 1)(s − 2)}

Problem 150. L(y) + L(t sin(t)) = F (s), F (s) = \frac{s + 1}{(s − 1)(s + 2)^{2}(s^{2} + 1)}

Problem 151. \frac{s}{s − 1}L(y) = \frac{6}{s^{3}} − \frac{2}{s + 1}

Special Methods. Solve for y(t) in the following equations using properties of Laplace transforms.

Problem 152. L(y) = L(tH(t − 1))\frac{s}{(s − 1)(s + 2)}

Problem 153. L(y′′(t)) + sy(0) + y′(0) = \frac{s}{s + 1}

Problem 154. \frac{d^{2}}{ds^{2}} L(y) − \frac{d}{ds}L(y) = F (s), F (s) = \frac{s − 1}{(s + 1)(s^{2} + 2s + 2)}

Problem 155. L(y) = L(t exp(−t))\frac{1}{s} arctan(s)

Convolutions and Systems. Solve for the unknowns in the following equations using Laplace transform theory.

Problem 156. L(∫_{0}^{t} sin(t − r)y(r)dr) = \frac{s}{(s^{2} + 1)^{2}}

Problem 157. L(x) + sL(y) = 1, (s − 1)L(x) − sL(y) = 0

Problem 158. L(∫_{0}^{t} y(r)dr) + 20L(y) = F (s), F (s) = \frac{1}{s^{2}} − \frac{exp(−s)}{s^{2}} − \frac{exp(−s)}{s}

Problem 159. L(x) + L(y) = \frac{s}{(s − 1)(s + 1)}, L(x) − L(y) = \frac{s}{(s − 2)(s + 2)}

Mutual Inductance. Two electric circuits are coupled inductively as shown in the figure below.

.......................................................................................................................  .............................................................................................................................................................................................................................................................................. ...................................................................................................................................................................................................................................  ................ ................  ............................................................................................................................................................................................................................................................................................. .............................................................................................................................................................................................................................................................................................  ...................................................................................................................................................................................................................................................................................................................... ........................................................................................................................................................................................................................................................................................................................  Secondary Primary E  R_{2} R_{1}  L_{2} I_{2} I_{1} M L_{1}

Mutual Inductance

Currents I_{1} and I_{2} are as shown. The number M is called the mutual inductance of the coupled inductors. Let L_{1} = L_{2} = 1, R_{1} = 1, R_{2} = 0.25, M = 0.5, E = 120. The model for the currents is:

L_{1} \frac{dI_{1}}{dt} + R_{1}I_{1} + M \frac{dI_{2}}{dt} = E,

L_{2} \frac{dI_{2}}{dt} + R_{2}I_{2} + M \frac{dI_{1}}{dt} = 0.

Problem 160. Show by a mathematical analysis that these two differential equations can be

re-written as

x′(t) = ax(t) + by(t) + f, y′(t) = cx(t) + dy(t) + g,

where x(t) = I_{1}(t), y(t) = I_{2}(t), K = L_{1}L_{2}−M^{2}, a = −R_{1}L_{2}/K, b = MR_{2}/K, f = EL_{2}/K, c = MR_{1}/K, d = −R_{2}L_{1}/K,g = −ME/K.

Problem 161. Solve the differential equations with zero initial data to obtain

x(t) = 120 + \frac{360}{\sqrt{13}} exp(−5t/6) sinh(\sqrt{13}t/6) − 120 exp(−5t/6) cosh(\sqrt{13}t/6),

y(t) = \frac{−480}{\sqrt{13}} exp(−5t/6) sinh(\sqrt{13}t/6).

Problem 162. How many seconds are required for the primary current to reach 99% of its steady state value of 120? How many seconds for the secondary current to reach 1% of its maximum excursion −39.23?

Transmission Line Voltage. In the LRC circuit of the figure below, the input voltage E(t) represents a signal, delivered to the circuit from a transmission line via a digital to analog converter. The goal is to compute the voltage across the capacitor

V (t) = q(t)/C

where q(t) is the charge on the capacitor of capacitance C.

.................. .................................................................................................................................................................................... ................................................................................................................................................................................................ ..................  .................. ................................................................................................................................................................  ................................................................................................................................................................................................................................................................................................................. q(t) C R E(t)  L V (t) i

Transmission Line

The governing differential equation for the circuit is (i = dq/dt):

L\frac{d^{2}q}{dt^{2}} + R\frac{dq}{dt} + \frac{q}{C} = E(t),

The assumed initial conditions on q are q(0) = q′(0) = 0.

Problem 163. Show that V (t) satisfies the differential equation and initial conditions

L\frac{d^{2}V}{dt^{2}} + R\frac{dV}{dt} + \frac{V}{C} = \frac{E(t)}{C}, V (0) = V′(0) = 0.

Problem 164. Let k(t) be the solution of the homogeneous differential equation determined by the Transfer Function of the circuit:

L(k(t)) = \frac{1}{Ls^{2} + Rs + 1/C} .

Assume L = 1, R = 2, C = 0.5. Find k.

Problem 165. Show by a transform theory analysis that V (t) is given by the convolution integral

V (t) = \frac{1}{C} ∫_{0}^{t} k(t − x)E(x)dx.

Problem 166. Select parameters L = 1, R = 2, C = 0.5 and let the emf E(t) = 10 exp(−t) sin(5t). Find V explicitly by integration of the convolution integral.

Problem 167. Numerical Methods for Mathematical Functions. The differential

equation to be solved is

y′ = y^{2} + g(x), y(0) = 0, 0 ≤ x ≤ 1.

Solve it for

g(x) = 2 cos(2 π x)π − cos(π x)π − 2 + cos^{2}(π x) + 2 sin(π x) sin(2 π x) + cos^{2}(2 π x)

by the Heun (Improved Euler) method and determine numerically the zeros and critical

points of y on [0, 1]. Approximate answers: zeros at 0, 0.3, 1 and critical points at 0.2, 0.7.

Numerical Methods for Experimental Functions.

Problem 168. Solve the differential equation

y′ = y^{2} + g(x), y(0) = 0, 0 ≤ x ≤ 1

numerically and determine the maximum and minimum values of y on [0, 1]. Assume the function g is given by

g(x) = (−24t + 3)H(t) + (21t − 7)H(t − \frac{1}{3}) + (51t − 34)H(t − \frac{2}{3})

This representation is the linear interpolant of the four data points

(0, 3), (1/3, −5), (2/3, −6), (1, 10).

The Heaviside function H is defined by

H(t − a) = \{ \begin{matrix}1 & t ≥ a \\ 0 & t < a \end{matrix}

‘ Problem 169. Solve the differential equation

y′ = y^{2} + g(x), y(0) = 0, 0 ≤ x ≤ 1

numerically and determine the maximum and minimum values of y on [0, 1]. Assume the function g is given by the linear interpolant of 11 data points (X_{0}, Y_{0}), …, (X_{10}, Y_{10}). The X-values start at 0.0 and stop at 1.0 in steps of 1/10 while the Y values are given in the table below:

    3.141592654,  2.017659785, -0.7319568225, -3.808368093, -6.185975224,
   -7.283185308, -6.480432161, -3.192887965,   2.115178899,  7.266781595,
    9.424777962

Definition. Assume given data points (X_{0}, Y_{0}), …, (X_{n}, Y_{n}) with X_{0} < ⋯ < X_{n}. Let H be Heaviside’s unit step function. Define the linear staircase function K_{1} by

K_{1}(t) = ∑_{k=1}^{n−1} H(t − X_{k})

Define the straight line function P_{1} and the linear interpolation function g of the data points {(X_{k}, Y_{k})}_{k=0}^{n} by

P_{1}(t, k) = \frac{Y_{k+1} − Y_{k}}{X_{k+1} − X_{k}} (t − X_{k}) + Y_{k}, g(t) = P_{1}(t, K_{1}(t))

Problem 170. (Numerical methods)Convert the differential equation u′′ + 3u′ + 2u = exp(t), u(0) = 1, u′(0) = 2 into the system of equations

x′ = y, y′ = −2x − 3y + exp(t), x(0) = 1, y(0) = 2

using the definition x = u, y = u′.

Problem 171. (Numerical methods)Solve the system

x′ = y, y′ = −2x − 3y + exp(t), x(0) = 1, y(0) = 2

by the Runge-Kutta method for x(1), y(1) using step size 0.2.

Problem 172. (Numerical methods)Solve the system of circuit equations

\frac{di_{1}}{dt} = −20i_{1} + 5i_{2} + 120, \frac{di_{2}}{dt} = 5i_{1} − 20i_{2}

with initial data i_{1}(0) = 0, i_{2}(0) = 1. Use Runge-Kutta and graph i_{1} and i_{2} on 0 ≤ t ≤ 2 in a single figure.

Problem 173. (Numerical methods)Solve by the Runge-Kutta method and graph on [0, π] the solution to the oscillator problem

u′′(t) + 4u(t) = g(t), u(0) = u′(0) = 0

where the external force g is given in terms of Heaviside’s function H by

g(t) = sin(t)H(t) + (cos(t) + 1 − sin(t))H(t − π/2)

Problem 174. (Numerical methods)Convert the Airy Differential Equation u′′ + tu = 0 with initial conditions u(0) = 0, u′(0) = 1 into a system of first order differential equations.

Problem 175. (Numerical methods)Convert the third order differential equation and initial

conditions

u′′′ − t^{2}u′′ + (u′)^{2}u = sin(t), u(0) = u′(0) = 0, u′′(0) = 1,

into a system of three first order differential equations.

Problem 176. (Numerical methods)Apply RK4 for systems to the Airy Differential Equation u′′ + tu = 0 with initial conditions u(0) = 0, u′(0) = 1 to produce the following table of values accurate to four digits:

\begin{matrix}x & y(x) & y′(x) \\ 0.0 & 0.0000000000 & 1.0000000000 \\ 0.4 & 0.3978699173 & 0.9787234992 \\ 0.8 & 0.7662804178 & 0.8329447612 \\ 1.2 & 1.0341745050 & 0.4643507191 \\ 1.6 & 1.1047663090 & −0.1469598482 \\ 2.0 & 0.8991799966 & −0.8834279315 \\ 2.4 & 0.4171370520 & −1.4736434800 \\ 2.8 & −0.2139584423 & −1.5744619860 \\ 3.0 & −0.5106490204 & −1.3612737990 \end{matrix}

Notes on the Problems

Maple. Below are a number of read commands which effectively enter source code into a maple program. These commands will be successful only in case the files exist.

Persons using maple on other hosts or on personal computers should copy required sources from math department directory /u/cl/maple. The copying is normally done with F T P. Personal computers can use a diskette copy made with mtools. Lab assistants in B129 can provide help

with ftp or mtools.

To use the source code below on personal computers or other hosts, change all references to

directory /u/cl/maple. For example, change read `/u/cl/maple/interpolation` to read
interpolation.

All the maple sources in /u/cl/maple have built-in help and examples. Use maple help to display the information.

Maple is supposed to have an initialization file present in your login directory. On the mathematics department computers, copy the maple initialization file to your root directory as follows:

       cp  -p   /u/cl/maple/sample.mapleinit   ~/.mapleinit

Problem Notes on 1 — 8.

The problem should be worked on paper without any help from maple. It is recommended that all problems be worked by hand in order to gain experience and intuition about differential equations. To check an answer using maple, apply the methods of the following example:

eq:= diff(x(t),t) = (1+t)/x(t);
ic:= x(0)=1;
dsolve({eq,ic},x(t));
#                                    d         1 + t
#                            eq := ---- x(t) = -----
#                                   dt          x(t)
#
#                                 ic := x(0) = 1
#
#                               2     1/2                 2     1/2
#              x(t) = - (2 t + t  + 1)   , x(t) = (2 t + t  + 1)
#

In some cases the integration required in the problem is difficult or even impossible. The services of maple can help in this situation, even though doing the problem by hand. For example, to integrate 1/(x^{2} + 3x + 2), apply the following maple command to compute the answer without the aid of an integral table (or partial fractions):

f:= x -> 1/(x^2+3*x+2);
int(f(x),x);
#                          - ln(x + 2) + ln(x + 1)

The answers:

Answer to 1. x(t) = \frac{ln((e^{2} − 1)e^{−2 t} + 1)}{2}

Answer to 2. \frac{x(t)^{3}}{3} + x(t) − \frac{t^{2}}{2} = 4/3 Answer to 3. y(x) − 5 ln(y(x) + 3) − x + 3 ln(x + 2) = 1 − 7 ln(2) Answer to 4. P (t) = e^{1−exp(−t)} Answer to 5. T (t) = 75 + 105 e^{− \frac{7 t}{10}}

Answer to 6. I(t) = −\frac{t^{4}}{1000} + \frac{t^{3}}{25} − \frac{3 t^{2}}{5} + 4 t

Answer to 7. y(x) = −\frac{4 − e^{ln( \frac{14}{11} )−3 x}}{1 − e^{ln( \frac{14}{11} )−3 x}} = \frac{44e^{3x} − 14}{14 − 11e^{3x}}

Answer to 8. y(x) = \frac{1}{2}\sqrt{8 e^{cos(x) sin(x)+x} − 4}

Some remarks are in order about problem 4, since certain maple errors are possible for ∫ dP/(P (1− ln P)). The answer is returned as − ln(1 − ln(P)), instead of − ln |1 − ln P |. If dsolve is used, then an incorrect term ln(ln P − 1) appears in the answer, which disallows P (0) = 1. It suffices to say that all troubles are avoidable by hand calculation.

Problem Notes on 9 — 16.

When solving linear equations y′ + p(x)y = q(x), apply the following working rule:

Let = ∫ p(x)dx.

Multiply the differential equation by exp( ).

The two-termed expression exp( )[y′ + p(x)y] is replaced by the single expression (exp( )y)′.

Example. Solve y′ + (1/x)y = 2. Solution: Define = ∫ (1/x)dx = ln x. Multiply the DE by exp( ). Replace the left hand side of the equation by (exp( )y)′. Integrate the resulting equation (exp( )y)′ = 2 exp( ) as in separation of variables to find the solution y:

exp( )[y′ + (1/x)y] = exp( ) · 2 (exp( )y)′ = exp( ) · 2 exp( )y = ∫ exp( ) · 2dx xy = ∫ x · 2dx xy = x^{2} + C

y = x + \frac{C}{x}

The example appears to be a collection of tricks because of the simplification exp( ) = exp(ln x) = x. The answers: Answer to 9. y(x) = e^{− arctan(x)} Answer to 10. x(t) = −t − 1 + 2 e^{t}

Answer to 11. x(t) = \frac{sin(t)}{2} − \frac{cos(t)}{2} + \frac{e^{−t}}{2} Answer to 12. x(t) = \frac{t^{2}}{3} + \frac{2}{3 t}

Answer to 13. U(t) = \frac{A}{B} − \frac{e^{−Bt}A}{B}

Answer to 14. v(t) = \frac{mg}{k} − e^{− \frac{kt}{m}} mgk^{−1}

Answer to 15. y(x) = \frac{3 x}{2} − \frac{3}{2 x}

Answer to 16. w(x) = x^{5}e^{x} − x^{4}e^{x}

Problem Notes on 17 – 20.

Given an exact differential equation M(x, y)dx + N(x, y)dy = 0, the implicit form of the general solution u(x, y) = c can be written as

∫_{x_{0}}^{x} M(t, y_{0})dt + ∫_{y_{0}}^{y} N(x, s)ds = c.

The point (x_{0}, y_{0}) in this formula is chosen in the domain of continuity of the functions M, N and c is a constant. If an initial condition is given, then c is zero and x_{0}, y_{0} are determined from the initial condition. Otherwise c is arbitrary and x_{0}, y_{0} are given specific values unrelated to initial conditions, usually x_{0} = y_{0} = 0 or x_{0} = y_{0} = 1, dictated by the domain of M and N.

Example. Given an exact equation Mdx + Ndy = 0 with initial condition y(1.1) = 4.6, then let c = 0, x_{0} = 1.1, y_{0} = 4.6 and use the formula

∫_{1.1}^{x} M(t, 4.6)dt + ∫_{4.6}^{y} N(x, s)ds = 0

Example. Evaluate ∫_{0}^{y} N(x, s)ds, given N = 2x + 5y.

Solution: Such integrals are called parametric integrals because the integrand and the limits of integration depend upon parameters x and y. The notation is succinct and historically troublesome.

If N(x, y) = 2x + 5y, then N(x, s) = 2x + 5s. The integral to be evaluated is ∫_{0}^{y}(2x + 5s)ds. The presence of the variable x in the integrand is unusual because in this context x is a constant! The integration is on variable s, therefore s is removed by the integration and ultimate substitution of s = 0 and s = y. The integral evaluates to 2xs + 5s^{2}/2 |_{s=0}^{s=y} = 2xy + 5y^{2}/2. In particular, the answer is a function of x and y and the variable s is not present due to integration.

Example. Solve the exact equation (x + y)dx + (x + 3y)dy = 0 with y(0) = 1.

Solution: The equation has the form Mdx + Ndy = 0 where M = x + y, N = x + 3y. It is exact because M_{y} = N_{x} = 1. Let x_{0} = 0, y_{0} = 1 and apply the formula:

∫_{0}^{x} M(t, 1)dt + ∫_{1}^{y} N(x, s)ds = 0

∫_{0}^{x}(t + 1)dt + ∫_{1}^{y}(x + 3s)ds = 0
(t^{2}/2 + t) |_{t=0}^{t=x} + (xs + 3s^{2}/2) |_{s=1}^{s=y} = 0
(x^{2}/2 + x) + x(y − 1) + 3y^{2}/2 − 3/2 = 0

The last equation is the implicit form of the solution.

The explicit form is obtained by solving this equation for y in terms of x. By the quadratic formula applied to Ay^{2} + By + C = 0 with A = 3/2, B = x, C = x^{2}/2 − 3/2 the explicit answer is

y(x) = \frac{−B ± \sqrt{B^{2} − 4AC}}{2A} = \frac{−x + \sqrt{−2x^{2} + 9}}{3}

The plus sign was selected for ± because of the requirement y(0) = 1.

If the explicit form cannot be obtained, then only the implicit form is reported.

The ability of maple to solve differential equations includes the class of exact equations. It is impossible to inform maple of the classification because it uses a fixed algorithm which cycles through the various classifications, looking for a solution. To solve an exact equation using maple, apply the ideas of this example:

eq:=(x+y)+ (x+3*y)*diff(y(x),x)=0:
ic:=y(0)=1:
dsolve({eq,ic},y(x));
#
#                                       1/2     2       1/2
#                 y(x) = - 1/3 x + 1/3 2    (- x  + 9/2)   ,
#
#                                           1/2     2       1/2
#                     y(x) = - 1/3 x - 1/3 2    (- x  + 9/2)
#
dsolve(eq,y(x));
#                            2                    2
#                       1/2 x  + x y(x) + 3/2 y(x)  = _C1
#

One defect in the maple answer is that the initial condition y(0) = 1 is ignored, causing two solutions to be presented. The presence of an initial condition sometimes prevents maple from finding a solution, therefore the second syntax used above is useful and in some cases necessary.

Answer to 17. \frac{5 x^{2}}{2} + 4 y(x)x − 2 y(x)^{4} = c

Answer to 18. cos(x) sin(y(x)) − ln(cos(x)) = c

Answer to 19. ln(x)e^{y(x)x} − y(x) = c

Answer to 20. x − 3 ln(x) + y(x)x + y(x) − 3 ln(y(x)) = c

Problem Notes on 21 – 28.

The possible answers to classification are:

separable and not linear linear and not separable linear and separable exact not separable, not linear, not exact

Separable and linear equations can be manipulated to an exact form in a standard way, therefore an exactness test is not performed for these equations. Answers like separable and exact or linear and exact are not given.

Test for separable first. Every separable equation dy/dx = F (x)G(y) can be put into exact form F (x)dx − (1/G(y))dy = 0. Test for linear second. Every linear equation y′ + p(x)y = q(x) can be multiplied by an integrating factor (x) = exp(∫ p(x)dx) to produce the exact equation [(p(x)y − q(x)] (x)dx + (x)dy = 0. Finally, if not separable and not linear, then test for exact. It is possible for an equation to be exact but not linear or separable.

An equation is separable if the y’s can be put on the left and the x’s on the right. Almost. This working rule seems to apply to y′ = x + exp(y) by writing it in the form y′ − exp(y) = x. More precisely, the working rule says that all y’s are left, all x’s are right and y′ is a factor on the left. Therefore, y′ = x + exp(y) is not separable.

To convert an equation like y′ = −x + exp(y) to the form Mdx + Ndy = 0, write y′ as dy/dx and then formally multiply by dx to clear fractions: (x − exp(y))dx + dy = 0. This equation is not separable, not linear and not exact (M_{y} ≠ N_{x}). An exactness test normally applies to the differential equation in its given form, without further manipulation.

The solutions returned by maple are as follows:

Answer to 21. x(t) = ln(exp(−t + ln(exp(1) − 1)) + 1)

Answer to 22. 2x^{5}(t) + 10x(t) − 5t^{2} − 10t − 12 = 0

Answer to 23. x(t) = −t − 1 + 2 exp(t)

Answer to 24. x(t) = −2 cos(t) + 2 sin(t) + 2 exp(−t)

Answer to 25. x(t) = (1/9)t^{5}

Answer to 26. y(x) = \frac{1}{1 + x^{2}}

Answer to 27. x exp(y(x)) − y(x) = 5

Answer to 28. y(t) = tanh(t) The maple commands that produced the above results are below. The trouble with problem 27 seems to be the initial data, because elimination of the initial data produces a usable answer. The answer in problem 22 is reported with a RootOf syntax, which has to be written again as an implicit equation.

eq:=diff(x(t),t) = -1 + \exp(-x(t)): ic:=x(0)=1:
dsolve({eq,ic},x(t));
eq:=diff(x(t),t) = (t+1)/((x(t))^4+1): ic:=x(0)=1:
dsolve({eq,ic},x(t));
eq:=diff(x(t),t) = x(t) + t: ic:=x(0)=1:

dsolve({eq,ic},x(t));
eq:=diff(x(t),t)  = -x(t) + 4*sin(t):  ic:=x(0)=0:
dsolve({eq,ic},x(t));
eq:=diff(x(t),t) = (-4/t)*x(t) + t^4: ic:=x(0)=0:
dsolve({eq,ic},x(t));
eq:=diff(y(x),x)*(1+x^2)+2*x*y(x): ic:=y(0)=1:
dsolve({eq,ic},y(x));
eq:=diff(y(x),x) *(x*exp(y(x))-1) + exp(y(x)) =0: ic:=y(5)=0:
dsolve({eq,ic},y(x)); dsolve(eq,y(x));
eq:=diff(y(t),t) = 1-y(t)^2: ic:=y(0)=0:
dsolve({eq,ic},y(t));

Problem Notes on 29 – 32.

The applied problems use terminology from circuit theory, biology and radioactive decay. Deep knowledge of these subjects is not required. Each application can be treated as a mathematics problem.

The circuit problems are first order linear equations to be solved by methods discussed previously. Answer to 29. Answer to 30. The exact answers are given by the maple solution:

eq:= 2*diff(x(t),t)+10*x(t)=200:
ic:=x(0)=0:
dsolve({eq,ic},x(t));
eq:= diff(x(t),t)+100*x(t)=sin(60*t):
ic:=x(0)=0:
dsolve({eq,ic},x(t),laplace);   # laplace often used for circuit equations
#
#                           x(t) = 20 - 20 exp(- 5 t)
#
#         x(t) = 3/680 exp(- 100 t) - 3/680 cos(60 t) + 1/136 sin(60 t)
#

The unusual use of laplace is necessary for this last example because of the strange treatment by maple V.2. If the laplace option is left off, then an answer is returned, but it appears to be excessively complicated and possibly even wrong. The symbolic answers needed for the last two problems can be obtained from maple as follows:

eq:= diff(x(t),t) = -1/3 + (1/12)*cos(Pi*(t-6)/24):
ic:=x(0)=600:
dsolve({eq,ic},x(t));
eq:= diff(x(t),t)=k*x(t):
ic:=x(0)=x0:
dsolve({eq,ic},x(t));
#                                     1/2                   1/2    1/2
#                     sin(1/24 Pi t) 2      cos(1/24 Pi t) 2      2
#    x(t) = - 1/3 t + ------------------- - ------------------- + ---- + 600
#                              Pi                    Pi            Pi
#
#                               x(t) = exp(k t) x0

Answer to 31. The problem requires that the solution be evaluated at t = 24: P (24) = 592.9gm.

Answer to 32. The symbolic answer contains two parameters k and x0. The value of k is determined from the half-life condition x(25) = 0.5x(0). The value of x0 is not determined nor is it necessary to know its value. The percentage remaining is calculated as 100 times the fraction remaining: 100x(32)/x(0). The final answer: 41.2% remains after 32 years.

Problem Notes on 33 – 36. A good reference for the Trapezoidal and Simpson rules is Kreyszig’s Advanced Engineering Mathematics, 7th edition, page 958 and page 961. Most modern calculus books have an elementary exposition on these two topics. It is expected that you can learn these ideas on your own from the references.

The exact answers are:

Answer to 33. y(x) = −x exp(−x) − exp(−x) + 1.

Answer to 34. y(x) = (1/2)(x^{2} − 1) ln(1 + x) − (1/4)x^{2} + (1/2)x + 1.

Answer to 35. y(x) = (−1/9)x^{9} + (1/6)x^{6}.

Answer to 36. (1/3) sin^{3}(πx).

The Trapezoidal and Simpson methods are already programmed into maple procedures in source file /u/cl/maple/numintegral. Copying the sources is recommended. The procedures are to be used as black boxes for numeric evaluation of an integral.

Example. Solve y′ = x exp(−x), y(0) = 0 exactly and also by the Trapezoidal Rule step 0.05. Graph the numerical solution.

Solution: The exact answer is found by integration using the Fundamental Theorem of Calculus: f(b) − f(a) = ∫_{a}^{b} f′(x)dx. In this application y(x) = ∫_{0}^{x} t exp(−t)dt = −x exp(−x) − exp(−x) + 1. Maple can be used as an integral table for such problems.

The numeric answer is found for 21 values and the plot is done in maple. The scheme used to generate the plot points is based upon the formula

∫_{a}^{b} f(x)dx = ∫_{a}^{c} f(x)dx + ∫_{c}^{b} f(x)dx

Here is how for the trapezoidal rule:

read `/u/cl/maple/numintegral`: # Get library routines trap() and simp()
f:= t -> t*exp(-t):             # Define function
a:=0.0: b:=1.0:                 # Define [a,b]
m:=20: h:=evalf((b-a)/m):       # Define step h
x0:=a: y0:=0.0:                 # Initial point
L:=[x0,y0]:                     # Make list of points for plotting
n:=10:                          # Use 10 trapezoidal divisions on [x0,x0+h]
for i from 1 to m do
  y0:=y0 + trap(x0,x0+h,n,f):   # Equals int(f,a..x0) + int(f,x0..x0+h);
  x0:=evalf(x0+h):              # (x0,y0)=next plot point
  L:=L,[x0,y0]:                 # Add new point to the list
od:
plot([ L ]);                    # Plot numeric solution (brackets required).

An explanation is required for the use of the symbol L above. In maple, the idea is to create a list, which is a special maple object consisting of strings separated by commas. In this case, the items of the list are themselves lists of two elements: [x, y]. Technically, L is a list of lists. The syntax L:=L,[x0,y0]: adds an item onto the current list. The plot package in maple can accept a list as input (i.e., an array of x and y values) for plotting. This particular plot option is very efficient and results in high speed plotting response.

A list like L above can be viewed by the maple command print(L). It is useful to plot the solution to a differential equation and then afterwards to examine closely the actual values used in the plot. This simplistic printing of the list does not align the columns or make a readable table, but often it is sufficient. The production of a printed table is made as follows. The table can be printed as it appears or saved in a file and edited in emacs to obtain the column alignment.

LL:=[ L ]:                      # Make an array from list L
N:=nops(LL):                    # Get size of list
for i from 1 to N do            # Print N rows in a table
print(LL[i][1],LL[i][2]);     # LL[i] is a list [x,y]
od:                             # LL[i][1] is x, LL[i][2] is y

The syntax is the same for Simpson’s Rule: replace trap by simp in the above example and change n to n = 5 (half the number of steps). It is essential to understand the maple code in the file numintegral. Study it in source form until you can understand each line.

# Trapezoidal rule for the integral of f(x) over [a,b] with n divisions.
# Here f:=proc(x) is a maple procedure for f(x), [a,b]
# is the interval, n is the number of subdivisions.
#
trap := proc(a,b,n,f)             # Trapezoidal Rule
local s,h,i:                      # Define temporary variables.
s:=evalf((f(b)+f(a))/2.0):        # s=accumulator
h:=(b-a)/n:                       # Define step size.
for i from 1 to n-1 do            # Add function values to the
  s:=evalf(s+f(evalf(a+i*h))):    #  accumulator for interior
od:                               #  divisions of [a,b].
s:=evalf(s*h):                    # Multiply accumulator by h.
RETURN(s):                        # Return the answer.
end:
#
# Simpson's rule for the integral of f(x) over [a,b] using 2*m
# divisions. Here f:=proc(x) is a maple procedure for f(x), [a,b]
# is the interval, m is half the number of subdivisions.
#
simp := proc(a,b,m,f)                     # Simpson's rule
local h,h2,s,i:
h:=(b-a)/(2*m): h2:=h+h: s:=0.0:          # Define step, accumulator.
for i from 1 to m do                      # Collect function values
  s:=evalf(s+f(evalf(a+i*h2-h))):         #  over odd interior
od:                                       #  divisions of [a,b].
s:=s*2.0:                                 # Double accumulator.
for i from 1 to m-1 do                    # Collect function values
  s:=evalf(s+f(evalf(a+i*h2))):           #  over even interior
od:                                       #  divisions of [a,b].
s:=evalf((s*2.0 + f(a) + f(b))*h/3.0):    # Double accumulator,  add x=a,
RETURN(s):                                #  x=b terms, multiply by h/3.
end:

Problem Notes on 37.

An analytical method is used to obtain what is often called a conservation law or an energy relation. Such a law leads directly to an explicit expression for y in terms of v.

The idea for the given problem is to let u(t) = 30 + 1.5x(t), then write the differential equation as (uu′)′ = −48u. Multiply this equation on both sides by uu′ and simplify to the form ww′ = −48u^{2}u′ where w = uu′.

Integrate ww′ = −48u^{2}u′ with respect to t over [0, T], then apply the power rule of integral calculus to obtain the conservation law. In the evaluations at t = 0 use the given information x(0) = 0, x′(0) = v. In terms of u and w, the evaluations at t = 0 are u(0) = 30 + 1.5x(0) = 30, u′(0) = 1.5x′(0) = 1.5v, w(0) = u(0)u′(0) = 45v. The evaluations at t = T are u(T) = 30 + 1.5x(T), w(T) = [30 + 1.5x(T)]1.5x′(T).

The maximum height y of the metal ball occurs at the time T when x′(T) = 0 and x(T) = y. Put these values into the conservation law to obtain the formula for y in terms of v.

The plot domain is [0, L] where L is between 6000 and 7000. To find the value of L, substitute y = 900 into the given equation and solve for v. The maple command solve can do this, but it is perhaps easier with a hand calculator.

A technique often used by engineers is to plot v versus y for larger and larger domains. When 900 appears on the ordinate axis, then the domain is large enough and an approximate value of L can be determined. This process can be very mysterious. A mathematical analysis is seen to be the only satisfactory method for determining the plot domain.

Answer to 37. The plot domain should be 0 to 6444. The function f to be plotted can be described by maple code

        g:= x-> (1/2)*(45*x)^2+16*(30)^3;
        f:= x -> (-30+(g(x)/16)^(1/3))/(3/2);

Problem Notes on 38.

The problem arose from the special kinetics problem of a rocket of weight 3000lb fired vertically upward, with assumed quadratic drag force coefficient of 0.008lb sec^{2}/ft^{2}. The number 32 was idealized from g = 32.2 ft/sec^{2}. In the original problem both the position x(t) and velocity v(t) = x′(t) are important. Here the velocity is investigated.

The use of maple on this problem is encouraged, especially to evaluate the integrals in the method of separation of variables, and finally, to check the answer using dsolve.

The differential equation v′ + (24/32)v^{2} = 32, v(0) = 0, has variables separable, and an explicit solution is possible for the velocity v(t).

Answer to 38. The answer for the velocity involves the tanh function defined by tanh(w) = (e^{w} − e^{−w})/(e^{w} + e^{−w}). Integration carried out by hand can produce other forms of the answer.

Problem Notes on 39 – 46.

The method for each problem is identical: (1) Substitute y = e^{λx} into the differential equation to obtain the characteristic equation in the variable λ, (2) Solve for the two roots λ, (3) apply the recipe (see page 75, Kreyszig’s Advanced Engineering Mathematics, 7th edition) to write down the general solution c_{1}y_{1} + c_{2}y_{2}. The answers:

Answer to 39. y_{1} = 1, y_{2} = x Answer to 40. y_{1} = 1, y_{2} = exp(−x) Answer to 41. y_{1} = cos(x), y_{2} = sin(x) Answer to 42. y_{1} = exp(−x), y_{2} = exp(−2x) Answer to 43. y_{1} = exp(−x), y_{2} = x exp(−x) Answer to 44. y_{1} = exp(−x/2) cos(\sqrt{3}x/2), y_{2} = exp(−x/2) sin(\sqrt{3}x/2) Answer to 45. y_{1} = cosh(3x), y_{2} = sinh(3x) Answer to 46. y_{1} = cos(x/3), y_{2} = sin(x/3) Problem Notes on 47 – 54. The method for each problem is identical: (1) Substitute y = x^{m} into the differential equation to obtain the Euler auxiliary equation in the variable m,

(2) Solve for the two roots m, (3) apply the recipe (see page 90, Kreyszig’s Advanced Engineering Mathematics, 7th edition) to write down the general solution c_{1}y_{1} + c_{2}y_{2}. The answers: Answer to 47. m(m − 1) = 0, y_{1} = 1, y_{2} = x. Answer to 48. m^{2} = 0, y_{1} = 1, y_{2} = log(x). Answer to 49. m^{2} − 2m = 0, y_{1} = 1, y_{2} = x^{2}. Answer to 50. m(m − 1) − 2m + 2 = 0, y_{1} = x, y_{2} = x^{2}. Answer to 51. 4m(m − 1) + 1 = 0, y_{1} = \sqrt{x}, y_{2} = \sqrt{x} log(x). Answer to 52. m(m − 1) − m + 2 = 0, y_{1} = x cos(log(x)), y_{2} = x sin(log(x)). Answer to 53. 4m(m − 1) − 4m + 3 = 0, y_{1} = \sqrt{x}, y_{2} = x\sqrt{x}. Answer to 54. m(m − 1) − 5m + 9 = 0, y_{1} = x^{3}, y_{2} = x^{3} log(x). Problem Notes on 55. Given y_{1} = x^{3}, y_{2} = x^{4} and the Euler differential equation x^{2}y′′ − 6xy′ + 12y = 0, to verify that y_{1}, y_{2} are linearly independent solutions on (0, ∞) it is required to (1) show y_{1} and y_{2} are solutions and (2) show they are linearly independent. The verification of (1) is by direct substitution of the given functions into the differential equation. The verification of (2) is based upon the theory of the Wronskian determinant: the functions are independent on (a, b) if they are continuously differentiable on (a, b) and their Wronskian determinant is nonvanishing for at least one value of x in (a, b). Problem Notes on 56. Given x^{3}, |x^{3}| and the Euler differential equation x^{2}y′′−4xy′+6y = 0, to verify that y_{1}, y_{2} are linearly independent solutions on (−∞, ∞) is a repetition of the method of problem 55. To prove the Wronskian determinant is identically zero requires a computation on (0, ∞) for x^{3} and x^{3} and a separate computation on (−∞, 0) for x^{3} and −x^{3}. At x = 0 the Wronskian determinant is defined and also zero. Does this violate the Wronskian Test? The test claims that dependent functions have zero Wronskian. It claims that the converse holds, if the coefficients a(x), b(x) of the differential equation y′′ + a(x)y′ + b(x)y = 0 are continuous on the given interval. Problem Notes on 57. Let x_{1}(t) be the response to external force g_{1}(t) in the pendulum equation x′′(t) + a^{2}x(t) = g_{1}(t), x(0) = x′(0) = 0. Suppose x_{2}(t) is the response for force g_{2}(t), x′′(t) + a^{2}x(t) = g_{2}(t), x(0) = c, x′(0) = d. To verify that the sum of the responses y(t) = x_{1}(t)+x_{2}(t) satisfies y′′+a^{2}y = g_{1}(t)+g_{2}(t) compute the left side y′′+a^{2}y as x′′_{1}+a^{2}x_{1}+x′′_{2}+a^{2}x_{2}. Use the differential equations for x_{1} and x_{2} to simplify to g_{1}+g_{2}, which is g. The initial conditions y(0) = c, y′(0) = d are treated similarly.

Problem Notes on 58. Given xy′′ + y′ = 0 and the solution y_{1}(x) = ln(x), a basis of solutions on (0, ∞) by the method of reduction of order is y_{1}, y_{2} where y_{2}(x) = u(x)y_{1}(x) and u is determined by the method: u(x) = ∫ exp(− ∫ P)/y_{1}^{2}(x)dx where P is the coefficient of y′ in normal form. In this case, P = 1/x, and u(x) = −(ln(x))^{−1}.

Problem Notes on 59. Given (1 − x^{2})y′′ − xy′ = 0 and the solution y_{1}(x) = 1, to determine a basis of solutions by the method of reduction of order use the formula y_{2}(x) = u(x)y_{1}(x) where u(x) = ∫ exp(− ∫ −x/(1 − x^{2})dx)/(1)^{2}dx.

Problem Notes on 60 – 65.

This collection of problems should be worked by hand with an assist by maple. The purpose of the exercise is to gain intuition about the methods and the form of the answers, not to find the answer. Use maple for integration problems, complicated substitutions and for checking intermediate answers.

The variation of parameters formula solves for a particular solution y_{p}(x) to the nonhomogeneous differential equation y′′ + p(x)y′ + q(x)y = f(x). The steps for constant coefficient equations:

1. Determine a basis y_{1}, y_{2} of solutions to the homogeneous differential equation

y′′ + p(x)y′ + q(x)y = 0.

2. Compute the Wronskian

W (t) = y_{1}(t)y′_{2}(t) − y_{2}(t)y′_{1}(t).

3. Compute the function k(x, t), known as Cauchy’s kernel, given by the formula

k(x, t) = \frac{−y_{1}(x)y_{2}(t) + y_{2}(x)y_{1}(t)}{W (t)}

4. Determine by direct integration a particular solution y_{p}(x) by the formula

y_{p}(x) = ∫_{0}^{x} k(x, t)f(t)dt

In Kreyszig, the formula for y_{p} is given as an indefinite integration, equivalent to the following integration formula:

y_{1}(x) ∫ \frac{−y_{2}(t)f(t)}{W (t)}dt + y_{2}(x) ∫ \frac{y_{1}(t)f(t)}{W (t)}dt.

One formula is a rearrangement of the other, they are not different, except that the one using definite integration contains extra terms not normally kept in the other. It may help in special cases to choose one form over the other, because the integration effort may be reduced.

Answers produced with indefinite integration will contain arbitrary constants which multiply the functions y_{1}(x) and y_{2}(x). These terms are usually absorbed into the general solution of the homogeneous equation, c_{1}y_{1} + c_{2}y_{2}. Answers obtained with definite integration will satisfy zero initial state conditions y_{p}(0) = y′_{p}(0) = 0. Terms in the solution which also appear in c_{1}y+1+c_{2}y_{2} may likewise be absorbed into the homogeneous solution.

Example. Find Cauchy’s kernel k for the equation y′′ + y = sin(2x) and solve for the general solution by the method of variation of parameters.

Solution: The homogeneous solution is given by the recipe for constant equations: c_{1}y_{1} + c_{2}y_{2} where y_{1} = cos(x), y_{2} = sin(x). By the formula for W (t) and k(x, t),

W (t) = cos(t) cos(t) − sin(t)(− sin(t)) = 1,

k(x, t) = \frac{− cos(x) sin(t) + sin(x) cos(t)}{W (t)} = sin(x − t),

because of the trig identity

sin(a − b) = sin(a) cos(b) − cos(a) sin(b).

To find the variation of parameters solution y_{p}(x) apply the formula:

y_{p}(x) = ∫_{0}^{x} k(x, t)f(t)dt.

Then

y_{p}(x) = ∫_{0}^{x} sin(x − t) sin(2t)dt

which becomes

y_{p}(x) = −(1/3) sin(2x) + (2/3) sin(x).

The general solution is therefore y = y_{p} + c_{1}y_{1} + c_{2}y_{2}, which becomes

y = (−(1/3) sin(2x) + (2/3) sin(x) + c_{1} cos(x) + c_{2} sin(x).

It is normal to combine (2/3) sin(x) and c_{2} sin(x) into a single term (2/3 + c_{2}) sin(x) and then rename 2/3 + c_{2} as c_{2}.

The same work can be done with the help of maple, especially to do the time-consuming integration phase.

de:=diff(k(x),x,x) + k(x)=0:  # DE for Cauchy's kernel k
ic:=k(0)=0,D(k)(0)=1:         # Cauchy kernel initial data
dsolve({de,ic},k(x),laplace); # Find it
#  k(x) = sin(x)              # Cauchy Kernel y''+ y=0
K:=unapply(rhs("),x):         # Isolate func in answer
f:=x -> sin(2*x):             # Define RHS of DE
int(K(x-t)*f(t),t=0..x):      # Integrate
simplify(");                  # Simplify

The answers for y_{p} by the method of variation of parameters are (and these are different than the method of undetermined coefficient answers, in general):

Answer to 60. \frac{14}{9} − \frac{x}{3} − \frac{2 xe^{2 x}}{3} − \frac{4 e^{2 x}}{9} − \frac{47 e^{−x}}{36} + \frac{7 e^{3 x}}{36} Answer to 61. \frac{e^{x}}{6} − \frac{e^{−x}}{6} Answer to 62. \frac{x}{25} + \frac{8}{625} + \frac{7 e^{4 x} sin(3 x)}{1875} − \frac{8 e^{4 x} cos(3 x)}{625} Answer to 63. \frac{cosh(2 x)x}{4} − \frac{e^{2 x}}{16} + \frac{sinh(2 x)x}{4} + \frac{e^{−2 x}}{16} = \frac{1}{4}xe^{2x} − \frac{1}{16}e^{2x} + \frac{1}{16}e^{−2x} Answer to 64. \frac{8}{3} + \frac{8 cos^{2}(x)}{3} − \frac{16 cos(x)}{3} Answer to 65. \frac{sin(x)}{3} − \frac{sin(x) cos(x)}{3} The answers obtained by the method of undetermined coefficients may be, and generally are, somewhat different, because the above answers for y_{p} must satisfy zero state conditions y(0) = y′(0) = 0. To convert the variation of parameters answer for problem 65 into the answer for the method of undetermined coefficients, write

c_{1} cos x + c_{2} sin x\frac{8}{3} + \frac{8 cos(x)^{2}}{3} − \frac{16 cos(x)}{3}

(c_{1} − \frac{16}{3}) cos x + c_{2} sin x\frac{8}{3} + \frac{8 cos^{2}(x)}{3}

Now the constant c_{1} − 16/3 can be called c_{1}, and the answer for undetermined coefficients is

c_{1} cos x + c_{2} sin x\frac{8}{3} + \frac{8 cos^{2}(x)}{3}

The form of y_{p} selected according to the method of undetermined coefficients: Answer to 60. y_{p} = c_{1} + c_{2}x + c_{3}e^{2x} + c_{4}xe^{2x} Answer to 61. y_{p} = c_{1}e^{x} Answer to 62. y_{p} = c_{1} + c_{2}x Answer to 63. y_{p} = c_{1}xe^{2x} Answer to 64. Use sin^{2}(x) = (1 − cos(2x))/2. y_{p} = c_{1} + c_{2} cos(2x) + c_{3} sin(2x) Answer to 65. Use sin(x) cos(x) = (0.5) sin(2x). y_{p} = c_{1} cos(2x) + c_{2} sin(2x) Example. Substitute a + bx + cx^{2} into the differential equation y′′ + y = 1 + x^{2} and determine the coefficients a, b, c. Solution: By hand it goes as follows:

2c + a + bx + cx^{2} = 1 + x^{2}
a + 2c = 1, b = 0, c = 1
a = −1, b = 0, c = 1

To do the same by machine is more work. The substitution step is possible and reasonably easy to do in practice:

f:= x -> a + b*x + c*x^2:
eq:= diff(f(x),x,x) + f(x) = 1 + x^2:

The last line printed by maple is the substitution of f into the differential equation. The equations a+2c = 1, b = 0, c = 1 are difficult to isolate by machine. Once isolated, the solution is convenient:

sys:={a+2*c=1,b=0,c=1}:
solve(sys,{a,b,c});

Example. Find the general solution of the differential equation y′′ + y = 1 + x^{2}.

Solution: The homogeneous solution is determined from the characteristic equation λ^{2} + 1 = 0 by the recipe: y_{c} = c_{1} cos(x) + c_{2} sin(x). The method of undetermined coefficients suggests a trial solution of the form y_{p} = a + bx + cx^{2}. By the previous example, the coefficients a, b, c are determined by substitution, giving the formula y_{p} = −1 + x^{2}. Therefore, the general solution y is determined by the superposition principle y = y_{c} + y_{p}:

y = c_{1} cos(x) + c_{2} sin(x) + x^{2} − 1.

The answer is checked using maple:

ODE:=diff(y(x),x,x)+y(x) = 1+x^2:
ANS:=y(x):
dsolve(ODE,ANS,laplace);

The answer in maple will involve y(0) and D(y)(0), which should be viewed as arbitrary constants.

Problem Notes on 66 – 71.

This collection of mathematical problems is designed to be solved from the theory with minimal knowledge of the applied material. The emphasis is on method, terminology and interpretation of the model. The tools used to solve the problems are variation of parameters and the method of undetermined coefficients. All equations have constant coefficients. The answers:

Answer to 66. The Cantilever Equation

EIy′′ = \frac{w(L − x)^{2}}{2} , y(0) = y′(0) = 0

has solution given by the variation of parameters formula because the initial conditions are zero state. The Cauchy kernel is (x − t)^{2}/2. The solution is

y(x) = \frac{wx^{4}}{24 EI} − \frac{wx^{3}L}{6 EI} + \frac{wx^{2}L^{2}}{4 EI}

Any work with this problem using maple should change E and I to lower case to avoid conflict with maple’s symbols E and I (2.7128 and \sqrt{−1}).

Answer to 67. The Overdamped Spring Equation

x′′(t) + 3x′(t) + 2x(t) = 0

is homogeneous and therefore the standard recipe applies to obtain the general solution c_{1}e^{−t} + c_{2}e^{−2t}.

Answer to 68. The Critically Damped Spring Equation

x′′(t) + 16x′(t) + 64x(t) = 0

can be solved by the recipe: c_{1}e^{−8t} + c_{2}te^{−8t}. The classification is changed if the constant 16 is replaced by 16.1 because the characteristic equation for 16.1 has distinct roots. Answer to 69. The Underdamped Spring Equation

x′′(t) + 2x′(t) + 2x(t) = 0

has solution

c_{1}e^{−t} cos(t) + c_{2}e^{−t} sin(t)

because the roots of the characteristic equation are −1 + i and −1 − i. Answer to 70. The forced, damped LRC circuit equation

(0.2)x′′ + (1.2)x′ + 2x = 4 sin(4t),
x(0) = x′(0) = 0

has by the method of undetermined coefficients the solution

x(t) = \frac{160 e^{−3 t} sin(t)}{51} + \frac{40 e^{−3 t} cos(t)}{51}

− \frac{10 sin(4 t)}{51} − \frac{40 cos(4 t)}{51}

The transient terms are those which tend to zero at infinity: they contain e^{−3t}. The remaining terms are the steady state terms. Steady state:

(−10/51) sin(4t) + (−40/51) cos(4t).

Transient:

(160/51)e^{−3t} sin(t) + (40/51)e^{−3t} cos(t).

The maple code for the full solution:

de:=(2/10)*diff(x(t),t,t)+
    (12/10)*diff(x(t),t)+2*x(t)=4*sin(4*t):
ic:=x(0)=0,D(x)(0)=0:
simplify(dsolve({de,ic},x(t),laplace));

Answer to 71. The forced, damped torsional motion equation Iθ′′ + cθ′ + kθ = T (t), θ(0) = 0, θ′(0) = 0 for I = 2, c = 2, k = 3, T (t) = sin(t) has by the method of undetermined coefficients the solution

x(t) = \frac{2 e^{− \frac{t}{2}} cos(\frac{\sqrt{5}t}{2})}{5} + \frac{sin(t)}{5} − \frac{2 cos(t)}{5}

If using maple on this problem, then replace I by i to avoid conflict with maple’s \sqrt{−1}. Problem Notes on 72. The RLC circuit equation

LQ′′ + RQ′ + \frac{1}{C} Q = e_{1} cos(wt) + e_{2} sin(wt),

Q(0) = q_{0}, Q′(0) = q_{1},

can be solved abstractly by considering the three cases for the discriminant of the characteristic equation. This is not suggested, because the algebra is especially complicated, making it difficult to gain any intuition about the problem. The particular cases of this problem cover the three possibilities for the characteristic equation and produce the three possible types of answers. The plots are interesting.

Each of the three data sets generates an example like the one below. The most challenging data set is the one in which the characteristic equation has complex roots.

Example. Solve

Q′′ + 3Q′ + Q = cos(10t) + sin(10t),
Q(0) = 0, Q′(0) = 1

exactly and plot the solution on [0, 5].

Solution: The characteristic equation is λ^{2}+3λ+1 = 0 with discriminant D = 5 and two distinct real roots (−3 + \sqrt{5})/2, (−3 − \sqrt{5})/2. The homogeneous equation, according to the recipe for constant equations, has general solution

y_{c} = c_{1}e^{(−3+\sqrt{5})t/2} + c_{2}e^{(−3−\sqrt{5})t/2}

Since the forcing term cos(10t) + sin(10t) corresponds to roots ±10i, which are not duplicated in roots of the characteristic equation λ^{2} + 3λ + 1 = 0, the method of undetermined coefficients applies to give the form of a trial particular solution as

y_{p} = a cos(10t) + b sin(10t)

with constants a and b to be determined. Substitution of the trial solution into the forced equation give the relation

−99a cos(10t) − 99b sin(10t) − 30a sin(10t) +30b cos(10t) = cos(10t) + sin(10t)

which in turn implies the algebraic equations

−99a + 30b = 1, 30a − 99b = 1

having solution a = −43/3567, b = −23/3567. At this point in the problem solution, the the general solution y = y_{c} + y_{p} is known with arbitrary constants c_{1} and c_{2} (as yet undetermined). The target solution is found by determining the constants c_{1} and c_{2} from the initial conditions y(0) = 0 and y′(0) = 1 (given as Q(0) = 0, Q′(0) = 1). The equations which result by substituting the general solution into the initial conditions are

c_{1}e^{0} + c_{2}e^{0} + a cos(0) + b sin(0) = 0
c_{1}e^{0}((−3 + \sqrt{5})/2) + c_{2}e^{0}((−3 − \sqrt{5})/2)
−10a sin(0) + 10b cos(0) = 1

or in reduced form

c_{1} + c_{2} = −a,

(−3 + \sqrt{5})/2)c_{1} + ((−3 − \sqrt{5})/2)c_{2} = 1 − 10b

The solution by Cramer’s rule is

c_{1} = \frac{43}{7134} + \frac{7723 \sqrt{5}}{35670}, c_{2} = \frac{\sqrt{5} (43 \sqrt{5} − 7723)}{35670}

and therefore the solution to the problem is approximately

Q(t) = 0.49017 exp(−0.38195t)−
0.47814 exp(−2.6181t)−
0.012055 cos(10t) − 0.0064480 sin(10t)

An account of the solution in maple appears below. The process is provided primarily for reference. It is not a suggestion for how to solve the problem, because of the complexity, but only as a means to check portions of the solution.

# Q'' + 3Q' + Q = cos(10t)+sin(10t), Q(0)=0, Q'(0)=1
a:='a': b:='b': c1:='c1': c2:='c2':                             # Variables
de:=diff(Q(x),x,x)+3*diff(Q(x),x)+Q(x)=cos(10*x)+sin(10*x):     # Forced differential equation
eval(subs(Q(x)=a*cos(10*x)+b*sin(10*x),de));                    # Stuff trial solution into DE
equ:={-99*a+30*b = 1,-30*a - 99*b = 1}:                         # Resulting equation
p:=solve(equ,{a,b});                                            # Solve for a, b
if lhs(p[1])='a' then                                           # Extract answers
  a:=rhs(p[1]): b:=rhs(p[2]):                                   # from the set of
else                                                            # answers given by
  a:=rhs(p[2]): b:=rhs(p[1]);                                   # maple's solve() function
fi:
f:= x -> a*cos(10*x) + b*sin(10*x):                             # Particular solution
g:= x -> c1*exp((-3+sqrt(5))*x/2)+c2*exp((-3-sqrt(5))*x/2):     # Homogeneous solution
eq1:=subs(x=0,f(x)+g(x)=0);                                     # Try to find c1, c2
eq2:=subs(x=0,diff(f(x),x)+diff(g(x),x)=1);                     # from Q(0)=0, Q'(0)=1
p:=solve({eq1,eq2},{c1,c2});                                    # Solve for c1, c2
if lhs(p[1])='c1' then                                          # Extract answers
  c1:=rhs(p[1]): c2:=rhs(p[2]):                                 # from the set of
else                                                            # answers given by
  c1:=rhs(p[2]): c2:=rhs(p[1]);                                 # maple's solve() function
fi:
evalf(g(x)+f(x),5);                                             # Print answer in readable form
plot(g(x)+f(x),x=0..5);                                         # Plot the solution on [0,5]

Problem Notes on 73.

(Convolution formula) The variation of parameters solution

x(t) = ∫_{0}^{t} k(t, r)g(r)dr

can be applied to x′′ + a^{2}x = g(t) because of the matching zero initial states.

The Cauchy kernel k is determined by solving the homogeneous problem

x′′ + a^{2}x = 0

for a basis of solutions x_{1}(t) = cos at, x_{2}(t) = sin at, then the Wronskian is

W (r) = a cos(ar) cos(ar) − a sin(ar)(− sin(ar))
= a

and the Cauchy kernel is

k(t, r) = \frac{cos(at) sin(ar) − sin(at) cos(ar)}{W (r)}

= \frac{1}{a} sin(at − ar).

This establishes the claimed formula. You are expected to fill in the general formulas that were used and add the missing details in this derivation. In particular, explain which trig identity was used in the last equality above, and the origin of the formula for k(t, r).

Problem Notes on 74.

(Resonance) The theory of mechanical resonance discusses the interaction between the frequency of the external force g(t) and the natural frequency of the unforced system x′′ + a^{2}x = 0. The term g(t) = 2a^{2} sin(at) has frequency a. The natural frequency of x′′ + a^{2}x = 0 is the natural frequency of its general solution c_{1} cos(at) + c_{2} sin(at), which is also a. Hence the frequencies are matched.

Matched frequencies cause the amplitudes of the solution x(t) to the forced system to grow without bound as t → ∞ because of the presence of unbounded terms t cos(at) or t sin(at) in the general solution.

The amplitude of x(t) = A sin(at)+B cos(at) is \sqrt{A^{2} + B^{2}}, therefore the solution x(t) = sin(at)− at cos(at) has amplitude \sqrt{1 + a^{2}t^{2}}. The latter has limit infinity at t = ∞.

To verify the claimed solution, substitute the proposed solution into the differential equation and initial conditions. To develop more insight into resonance, solve the problem from the point of view of variation of parameters. In this setting it is possible to see why a term t cos(at) must appear in the answer.

Problem Notes on 75.

(Wine glass experiment)

Suppose x(t) approximates the deflection from equilibrium of the rim of a wine glass subjected to sound waves from a loudspeaker. The Hooke’s law term ax(t) would be derived from elastic properties of the wine glass. The forcing term g(t) would be derived from the sound waves emitted by the loudspeaker. Explain in a short paragraph how the theory of resonance applies to design an experiment to break the wine glass.

The classical experiment was broadcast over public television in the Annenburg Project series on resonance. The film is available from video vendors and rental is possible.

In this series a loudspeaker is placed near a wine glass. The input to the loudspeaker is a sound wave of fixed frequency, essentially g(t) = A sin(bt) for some b. The value of b is adjusted by tuning a frequency generator and the value of A is adjusted from an audio amplifier. The natural frequency of the wine glass is a number a which is a property of the material. The height of the glass and the length of the stem also affect the vibrations of the glass.

When b is very close to a the vibrations of the rim of the wine glass grow until the elastic limits of the glass are exceeded. There is always damping present in such a system, however, it is negligible compared to the excursions of the glass rim caused by matching of the external force to the glass frequency. Indeed, each time the rim moves the external force (the sound wave) is there to assist the excursion, making larger and larger amplitudes.

Problem Notes on 76 – 79. The technique used to solve higher order constant equations depends entirely on being able to factor the characteristic equation. Cubics always have one real root, therefore long division by the corresponding factor produces a quadratic polynomial whose roots are the other two roots of the cubic.

Example. Solve λ^{3} + 3λ^{2} + 4λ + 2 = 0 for x.

According to the rational root theorem, the possible rational roots of the equation are 1, −1, 2 and −2, all of which are factors of the constant term (2) divided by factors of the leading coefficient (1 is the coefficient of λ^{3}).

By the signs of the terms, it follows that only negative roots are possible. So we try −1 and −2, finding out that −1 is a root.

Therefore, λ − (−1) is a factor. Dividing by λ + 1, using long division, gives quotient λ^{2} + 2λ + 2, which has two complex conjugate roots −1 ± i.

The complete set of roots is −1, −1 + i and −1 − i.

The answers can be checked in maple by using the solve command, as follows:

     solve(u^3+3*u^2+4*u+2=0,u);

The answers for the roots duplicate what was said above: −1, −1 + i and −1 − i.

Example. Solve y′′′ + 3y′′ + 4y′ + 2y = 0.

The basis elements of the general solution are generated from the roots. By the previous example, the roots are −1, −1 + i and −1 − i. The first root −1 gives e^{−x}. The complex conjugate roots give e^{−x} cos x and e^{−x} sin x. Therefore, the general solution is

y = c_{1}e^{−x} + c_{2}e^{−x} cos x + c_{3}e^{−x} sin x.

Answers can be checked in maple by using the dsolve command, as follows:

de:=diff(y(x),x,x,x)+
    3*diff(y(x),x,x)+4*diff(y(x),x)+
    2*y(x)=0:
dsolve(de,y(x),laplace);

The answer for the differential equation is expressed in terms of y(0), D(y)(0) and D(D(y))(0), which should be viewed as arbitrary constants.

Example. Solve λ^{3} + 3λ^{2} + 4λ + 1 = 0 for λ.

This example cannot be done by college algebra methods, because the real root is irrational. An approximation to the real root is obtained by Newton’s method as −0.31764. Applying long division to the cubic λ^{3} + 3λ^{2} + 4λ + 2 with divisor λ + 0.31764 gives the quadratic

λ^{2} + 2.68236λ + 3.14797517

which has complex roots approximately −1.341 + 1.162i and −1.341 − 1.162i. The complete set of roots is approximately −0.3176, −1.341 + 1.162i and −1.341 − 1.162i. They can be checked by maple as follows:

evalf(solve(u^3+3*u^2+4*u+1=0,u),5);

Example. Solve x^{4} − 10x^{3} + 35x^{2} − 50x + 24 = 0. The factors of 24 give the possible rational roots, and substitution reveals that 1, 2, 3 and 4 are roots. Fourth order equations have exactly 4 roots, therefore all roots have been found. Maple gives us an exact solution 1, 2, 3, 4 from the command

solve(x^4-10*x^3+35*x^2-50*x+24=0,x);

Example. Solve x^{5} + x^{4} − 7x^{3} − 3x + 2 = 0. The technique of finding rational roots suggests only the values ±1 and ±2, none of which are roots. A practical way to proceed is by invoking maple as follows:

allvalues(solve(x^5+x^4-7*x^3-3*x+2=0,x));

−3.255195196,
0.4605880040,
2.267993844,
−.2366933264 − .7294794757 i,
−.2366933264 + .7294794757 i.

Answer to 76. x(t) = C_{1} e^{−t} + C_{2} e^{−t} cos(t) + C_{3} e^{−t} sin(t) Answer to 77. x(t) = C_{1} e^{1.34t} + C_{2} e^{t} cos(\sqrt{3}t) + C_{3} e^{t} sin(\sqrt{3}t) Answer to 78. x(t) ≈ C_{1} e^{1.336788104t}+C_{2} e^{1.001605948t} cos(1.731737767t)+C_{3} e^{1.001605948t} sin(1.731737767t). Answer to 79. x(t) = C_{1} e^{7t/3} + C_{2} e^{\sqrt{3}t} + C_{3} te^{7t/3} + C_{4} e^{−\sqrt{3}t} Problem Notes on 80 – 83. Example. Solve y′′′ + 3y′′ + 4y′ + 2y = 0 and plot the solution with initial conditions y(0) = 1, y′(0) = 0, y′′(0) = 0. By a previous example the general solution is given by the formula

y = c_{1}e^{−x} + c_{2}e^{−x} cos x + c_{3}e^{−x} sin x.

We have to determine c_{1}, c_{2} and c_{3} satisfying the initial conditions y(0) = 1, y′(0) = 0, y′′(0) = 0. The equations for the constants are

c_{1} + c_{2}  = 1,
−c_{1} − c_{2} + c_{3} = 0,
c_{1}  − 2c_{3} = 0.

These equations can be solved with Cramer’s rule as follows:

D = | \begin{matrix}1 & 1 & 0 \\ −1 & −1 & 1 \\ 1 & 0 & −2 \end{matrix} | = 1,

D_{1} = | \begin{matrix}1 & 1 & 0 \\ 0 & −1 & 1 \\ 0 & 0 & −2 \end{matrix} | = 2,

D_{2} = | \begin{matrix}1 & 1 & 0 \\ −1 & 0 & 1 \\ 1 & 0 & −2 \end{matrix} | = −1,

D_{3} = | \begin{matrix}1 & 1 & 1 \\ −1 & −1 & 0 \\ 1 & 0 & 0 \end{matrix} | = 1,

c_{1} = \frac{D_{1}}{D} = 2, c_{2} = \frac{D_{2}}{D} = −1, c_{3} = \frac{D_{3}}{D} = 1.

Therefore, back-substitution of these values in the general solution gives the explicit solution

y = 2e^{−x} − e^{−x} cos x + e^{−x} sin x.

The plot command in maple is

f:=x->2*exp(-x) - exp(-x)*cos(x )
      + exp(-x)*sin(x);
plot(f(x),x=0..10);

Example. Plot the solution of y′′′ + 3y′′ + 4y′ + y = 0 with initial conditions y(0) = 1, y′(0) = 0, y′′(0) = 0.

The solution with the given initial conditions is obtained by specializing the values of c_{1}, c_{2} and c_{3} in the general solution (see the previous example)

\begin{matrix}y = c_{1}e^{−0.3176x} + c_{2}e^{−1.341x} cos(1.162x) \\ + c_{3}e^{−1.341x} sin(1.162x). \end{matrix}

The equations for the constants are

\begin{aligned}c_{1} & + & c_{2} & + & 0c_{3} & = & 1 \\ −0.3176c_{1} & + & (−1.341)c_{2} & + & (1.162)c_{3} & = & 0 \\ (0.3176)^{2}c_{1} & + & 0.448037c_{2} & − & 3.116484c_{3} & = & 0 \end{aligned}

They can be solved by Cramer’s rule to give approximate answers

\begin{aligned}c_{1} & = & 1.313203238, \\ c_{2} & = & −0.3132032378, \\ c_{3} & = & −0.00252340252. \end{aligned}

The function to be plotted is

\begin{matrix}g(x) = 1.3132e^{−0.3176x} − 0.31320e^{−1.341x} cos(1.162x) \\ − 0.0025234e^{−1.341x} sin(1.162x). \end{matrix}

Maple completes the plot problem by the command

g:= x->1.3132*exp(-0.3176*x)
       -0.3132*exp(-1.341*x)*cos(1.162*x)
       -0.0025234*exp(-1.341*x)*sin(1.162*x);
plot(g(x),x=0..10);

Much of the tedious part of this problem can be done with maple assist. The details:

f:=x->c1*exp(-0.3176*x)
      +c2*exp(-1.341*x)*cos(1.162*x)
      +c3*exp(-1.341*x)*sin(1.162*x);
sys:={f(0)=1,D(f)(0)=0,D(D(f))(0)=0};
solve(sys,{c1,c2,c3})  # Solve the system.
evalf(",5);            # Cut off at 5 digits.
subs(",f(x));          # Form evaluated equation.
g:=unapply(",x);       # Call it g(x).
plot(g(x),x=0..10);    # Plot the solution.

Example. Solve x′′′ − 4x′′ − 80x′ − 192x = 0, x(0) = 1, x′(0) = 0, x′′(0) = 0 using maple. Proceed as follows:

ODE:=(D@@3)(x)(t)-4*(D@@2)(x)(t)
     -80*D(x)(t)-192*x(t)=0;
IC:=x(0)=1,D(x)(0)=0,(D@@2)(x)(0)=0;
dsolve({ODE,IC},x(t));

x( t) = \frac{15}{16} e^{−4t} + \frac{1}{16} e^{12t)} + 3 e^{−4t} t

Answer to 80. x(t) = \frac{15}{16}e^{−4t} + \frac{1}{16}e^{12t} + 3te^{−4t} Answer to 81. x(t) = −\frac{1}{12}e^{t} + \frac{5\sqrt{273}}{156}e^{− \frac{t}{8}} sinh(\frac{\sqrt{273}t}{8}) + \frac{1}{12}e^{− \frac{t}{8}} cosh(\frac{\sqrt{273}t}{8}) Answer to 82. x(t) = 3e^{−3t} − 8e^{−2t} + 6e^{−t} Answer to 83. x(t) = te^{−t} + \frac{3}{2}e^{−t} − \frac{1}{2} cos t + \frac{1}{2} sin t

Problem Notes on 84 – 87.

Example. Solve the equation x′′ + 2t x′ + 2x = 0. The general solution of this equation is x = A_{1}x_{1} +A_{2}x_{2} for arbitrary constants A_{1} and A_{2}, where the functions x_{1} and x_{2} are linearly independent solutions obtained as series expansions by the power series method. The solution x_{1} will be shown below to be given by the formula

x_{1}(t) = ∑_{k=0}^{∞} \frac{(−t^{2})^{k}}{k!}

which by series formulas from calculus equals e^{−t^{2}}. The second solution x_{2}, given by

x_{2}(t) = ∑_{k=0}^{∞} \frac{k!(−4)^{k}}{(2k + 1)!} t^{2k},

is obtained similarly.

The trial solution will be x(t) = ∑_{k=0}^{∞} c_{k}t^{k}. The method: find formulas for x′ and x′′, then insert the formulas in the differential equation, obtain a recursion relation, solve the recursion and present the answer. Some details:

x(t) = ∑_{k=0}^{∞} c_{k}t^{k}  Trial solution.

x′(t) = ∑_{k=0}^{∞} kc_{k}t^{k−1} Differentiation Rule.

tx′(t) = ∑_{k=0}^{∞} kc_{k}t^{k}  Power Rule.

x′′(t) = ∑_{k=2}^{∞} k(k − 1)c_{k}t^{k−2} Differentiation Rule.

= ∑_{n=0}^{∞}(n + 2)(n + 1)c_{n+2}t^{n} Indexing rule with n = k − 2.

{∑_{n=0}^{∞}(n + 2)(n + 1)c_{n+2}t^{n} + \atop 2 ∑_{k=0}^{∞} kc_{k}t^{k} + 2 ∑_{k=0}^{∞} c_{k}t^{k} = 0}

Insert into x′′ +2tx′ +2x = 0. Use the sum rule.

∑_{k=0}^{∞} [(k + 2)(k + 1)c_{k+2} + 2kc_{k} + 2c_{k}] t^{k} = 0

(k+2)(k+1)c_{k+2}+2(k+1)c_{k} = 0 Vanishing Rule, k = 0, 1, 2, . . ..

c_{k+2} = (\frac{−2}{k + 2}) c_{k} Recursion. k = 0, 1, 2, . . ..

Basis solution x_{1}(t)

To solve the recursion relations and obtain x_{1}, set c_{0} = 1 and c_{1} = 0 (position and velocity).

Then:

c_{2n+1} = 0, n = 0, 1, 2, … By the recursion and c_{1} = 0.

c_{2} = −1 Recursion with k = 0, using c_{0} = 1.

c_{4} = \frac{−c_{2}}{2} = \frac{1}{2}

Recursion with k = 2, using c_{2} = −1.

c_{6} = −1/3!

Recursion with k = 4, using c_{4} = 1/2.

c_{2n} = (−1)^{n}/n!

General pattern. n = 0, 1, 2, …

The solution x_{1} is obtained by splitting the formal series into even and odd terms:

x_{1}(t) = ∑_{k=0}^{∞} c_{k}t^{k}  Trial solution.

{= ∑_{n=0}^{∞} c_{2n}t^{2n} \atop + ∑_{n=0}^{∞} c_{2n+1}t^{2n+1}} Even and odd terms.

= ∑_{n=0}^{∞} c_{2n}t^{2n}  {Because \atop c_{2n+1} = 0.}

= ∑_{n=0}^{∞} \frac{(−1)^{n}}{n!} t^{2n}  {Because \atop c_{2n} = (−1)^{n}/n!.}

Basis solution x_{2}(t)

The second solution x_{2}(t) is constructed with initial values c_{0} = 0 and c_{1} = 1. The analysis is basically unchanged, giving even coefficients c_{2k} = 0 for k = 0, 1, 2, … and odd coefficients c_{3} = −2/3, c_{5} = 2^{2}/3 · 5, c_{7} = −2^{3}/3 · 5 · 7. The general formula for the odd coefficients is

c_{2k+1} = \frac{(−4)^{k}k!}{(2k + 1)!}, k = 0, 1, 2, …

which gives the second basis solution

x_{2}(t) = ∑_{n=0}^{∞} \frac{n!(−4)^{n}}{(2n + 1)!} t^{2n+1}.

The recursion formulas and the truncated series x(t) (to 8 terms) are given below. Graphs of the solution generally use the formula for the truncated series.

Answer to 84. c_{n} = −\frac{c_{n−3}}{n(n − 1)},

c_{2} = 0, x(t) = t − \frac{1}{12}t^{4} + \frac{1}{504}t^{7}

Answer to 85. c_{n} = −\frac{(2n − 5)c_{n−2}}{n(n − 1)} ,

x(t) = t − \frac{1}{6}t^{3} + \frac{1}{24}t^{5} − \frac{1}{112}t^{7}

Answer to 86. c_{n} = − ∑_{j=0}^{n−2} \frac{c_{n−2−j}}{j!n(n − 1)},

x(t) = 1 − \frac{t^{2}}{2} − \frac{t^{3}}{6} + \frac{t^{5}}{40} + \frac{t^{6}}{80} + \frac{t^{7}}{315}

Answer to 87. c_{n} = −\frac{c_{n−4}}{n(n − 1)},

c_{2} = c_{3} = 0

x(t) = t − \frac{1}{20}t^{5}

Example. Check the series answer for x′′ + tx = 0, x(0) = 0, x′(0) = 1 using maple.

The maple commands:

Order:=8:
ODE:=diff(x(t),t,t)+t*x(t)=0:
IC:=x(0)=0,D(x)(0)=1:
ANS:=x(t):
dsolve({ODE,IC},ANS,series);

The answer generated is

x(t) = t − \frac{1}{12}t^{4} + \frac{1}{504}t^{7} + O (t^{8}) .

The notation O (t^{8}) stands for the sum of the series from the term containing t^{8} to infinity. It is read big Oh of t^{8}.

Problem Notes on 88. By the theory, trial solutions for P_{0} and P_{1} are respectively y = c and y = ax + b. Substitute these trial solutions into the differential equation and compute a, b, c. Then use P (1) = 1 to evaluate any unknowns.

Problem Notes on 89. The polynomials P_{0}(x) = 1 and P_{1}(x) = x are obtained from problem 88. Apply the Legendre polynomial recurrence relation to obtain P_{2}(x) to P_{8}(x).

Answer to 89. The first eight Legendre Polynomials are

P_{1}(x) = x

P_{2}(x) = \frac{3x^{2}}{2} − 1/2

P_{3}(x) = \frac{5x^{3}}{2} − \frac{3x}{2}

P_{4}(x) = \frac{35x^{4}}{8} − \frac{15x^{2}}{4} + 3/8

P_{5}(x) = \frac{63x^{5}}{8} − \frac{35x^{3}}{4} + \frac{15x}{8}

P_{6}(x) = \frac{231x^{6}}{16} − \frac{315x^{4}}{16} + \frac{105x^{2}}{16} − \frac{5}{16}

P_{7}(x) = \frac{429x^{7}}{16} − \frac{693x^{5}}{16} + \frac{315x^{3}}{16} − \frac{35x}{16}

P_{8}(x) = \frac{6435x^{8}}{128} − \frac{3003x^{6}}{32} + \frac{3465x^{4}}{64} − \frac{315x^{2}}{32} + \frac{35}{128}

Problem Notes on 90. The Legendre recurrence with x = 1 becomes (n + 1)u_{n+1} = (2n + 1)xu_{n} − nu_{n−1} where u_{n} = P (1). Given u_{0} = u_{1} = 1, then by solving this recursion successively, it follows that u_{n} = 1 for all n. The proper method for doing the proof is mathematical induction.

Problem Notes on 91. The graphs are produced by the maple command plot(f(x),x=-1..1), where f is one of the polynomial answers from problem 89.

Problem Notes on 92 – 101.

The main discussion about method appears below. The challenge is to learn the theory and fundamental methods. Answers can be obtained from maple by modifying the following example. A complete list of answers with additional notes appears below.

Example. Compute L(3te^{−t} sin(t) + cosh(t)).

Solution: The solution by hand is as follows:

L(3te^{−t} sin(t) + cosh(t)) – Original problem. L(3te^{−t} sin(t)) + L((e^{t} + e^{−t})/2) – Sum rule. 3L(te^{−t} sin(t)) + (1/2)L(e^{t} + e^{−t}) – Constant rule. 3L(te^{−t} sin(t)) + (1/2)L(e^{t}) + (1/2)L(e^{−t}) – Sum rule.

3L(te^{−t} sin(t)) + \frac{1}{2(s − 1)} + \frac{1}{2(s + 1)} – Exponential rule.

−3(d/ds)L(e^{−t} sin(t)) + \frac{1}{2(s − 1)} + \frac{1}{2(s + 1)} – Differentiation rule.

−3(d/ds) L(sin(t))|_{s→(s+1)} + \frac{1}{2(s − 1)} + \frac{1}{2(s + 1)} – Shifting rule.

−3(d/ds)\frac{1}{(s + 1)^{2} + 1} + \frac{1}{2(s − 1)} + \frac{1}{2(s + 1)} – Tables.

3\frac{2s + 2}{((s + 1)^{2} + 1)^{2}} + \frac{1}{2(s − 1)} + \frac{1}{2(s + 1)} – Derivative in s.

This answer can be left as it stands, or further algebra can be applied to make it a single fraction. The fractional answer appears below, as it emerges from maple. In maple the following code is used:

readlib(laplace):               # Load library
L:=laplace:                     # Abbreviation for laplace operator
L(3*t*exp(-t)*sin(t) + cosh(t),t,s);
#

The answer from maple is

\frac{14 s^{3} − 2 s + 14 s^{2} − 6 + s^{5} + 4 s^{4}}{(s^{2} + 2 s + 2)^{2} (s^{2} − 1)}

Answer to 92. L(3t + 2) = 3L(t) + 2L(1), linearity of the transform. Answer: \frac{2}{s} + \frac{3}{s^{2}}

Answer to 93. L(t^{2}) = L((−t)(−t)·1) = (d/ds)^{2}L(1), differentiation rule. Answer: \frac{1.1}{s} − \frac{4}{s^{2}} + \frac{2}{s^{3}}

Answer to 94. (t − 1)(t + 1)(t − 2) = (t^{2} − 1)(t − 2) = t^{3} − 2t^{2} − t + 2. Answer: \frac{2}{s} − s^{−2} − \frac{4}{s^{3}} + \frac{6}{s^{4}}

Answer to 95. L(5e^{1−t}) = 5L(e^{1}e^{−t}) = 5eL(e^{−t}). Answer: \frac{5 e^{1}}{s + 1}

Answer to 96. sinh(−t/3) = (1/2)(e^{−t/3} − e^{−(−t/3)}). Answer: \frac{9 s − 6}{9 s^{2} − 1}

Answer to 97. L((t^{2}−t)e^{−t}) = L(t^{2}−t) |_{s→(s+1)} , shifting theorem. Answer: \frac{2}{(s + 1)^{3}} −(s + 1)^{−2}

Answer to 98. L(e^{2t} cos(4t)) = L(cos(4t)) |_{s→(s−2)} , shifting theorem. Answer:

−\frac{10 s − 5 s^{2} − 24 + s^{3}}{(s^{2} − 2 s + 2) (s^{2} − 4 s + 20)}

Answer to 99. L(t^{2} cos(3t)) = (d/ds)^{2}L(cos(3t)), differentiation theorem. Answer:

\frac{8 s^{3}}{(s^{2} + 9)^{3}} − \frac{6 s}{(s^{2} + 9)^{2}}

Answer to 100. L(te^{−t} sin(2t)) = (−d/ds)L(e^{−t} sin(2t)) = (−d/ds)L(sin(2t)) |_{s→(s+1)} , differen-

tiation theorem, shifting theorem. Answer: \frac{4 s + 4}{((s + 1)^{2} + 4)^{2}}

Answer to 101. L(2te^{−3t} cosh(t)) = L((−1)(−t)(e^{−2t} + e^{−4t})) = (−d/ds)L(e^{−2t} + e^{−4t}), differ-

entiation theorem. Answer: \begin{matrix}(2 s + 6) & 2 \\ ((s + 3)^{2} − 1) & 2 \end{matrix} − \frac{2}{(s + 3)^{2} − 1}

Problem Notes on 102 – 109.

The answers and notes on method appear below.

Answer to 102. sin^{2}(3t) = 0.5 − 0.5 cos(6t) by a trigonometric identity. Answer: \frac{18}{s (s^{2} + 36)}.

The maple commands:

readlib(laplace): L:=laplace:
L((sin(3*t))^2,t,s);

Answer to 103. Let f(t) = sin(t)/t. Then L(tf(t)) = L(sin(t)) and by the differentiation rule (d/ds)L(f(t)) = −1/(1 + s^{2}). The function f is piecewise continuous and of exponential order, therefore by a theorem lim_{s→∞} L(f) = 0. Integrate over [t, ∞) to obtain an equation for L(f(t)). Answer: \frac{π}{2} − arctan(s)

Answer to 104. (1 − cosh(2t))/t can be transformed from the abstract result ∫_{t}^{∞} L(f)ds =

L(f(t)/t) (see above). Answer: \frac{ln(s − 2)}{2} − ln(s) + \frac{ln(s + 2)}{2} or ln \sqrt{1 − (4/s)^{2}}

Answer to 105. e^{t} + cos(t − π)H(t − π) is transformed by the standard theorem on Heaviside

functions: L(f(t − a)H(t − a)) = e^{−as}L(f(t)). Answer: (s − 1)^{−1} + \frac{e^{−sπ}s}{s^{2} + 1}

Answer to 106. H(t) + H(t − 1) + H(t − 2) + H(t − 3) is transformed by the theorem on

Heaviside’s function H: L(H(t − a)) = e^{−as}/s. Answer: \frac{1}{s} + \frac{e^{−s}}{s} + \frac{e^{−2 s}}{s} + \frac{e^{−3 s}}{s}

Answer to 107. | sin(2t)| is a periodic function of period T = π/2. Apply the theorem L(f) =

∫_{0}^{T} e^{−st}f(t)dt/(1 − e^{−sT}). Answer: \frac{2}{s^{2} + 4} \frac{1 + e^{−πs/2}}{1 − e^{−πs/2}}

Answer to 108. H(t) + H(t − 1) on [0, 2] is transformed in the same way as in problem 107 above using f(t) = H(t) + H(t − 1), T = 2. The evaluation of the integral is unusual: ∫_{0}^{2} H(t)e^{−st}dt = ∫_{0}^{2} e^{−st}dt because H(t) = 1 for t ≥ 0. Similarly, ∫_{0}^{2} H(t − 1)e^{−st}dt = ∫_{1}^{2} e^{−st}dt

because H(t − 1) = 0 for 0 ≤ t < 1. Answer: \frac{2 e^{−2 s} − 1 − e^{−s}}{(e^{−2 s} − 1) s}. The maple command:

(1/(1 - exp(-2*s))) *(int(exp(-s*x), x=0..2) +int(exp(-s*x), x=1..2));

Answer to 109. tH(t) + (1 − 2t)H(t − 1) on [0, 2] is similar to problem 108 above. Answer:

−\frac{e^{−2 s}s + e^{−2 s} + 1 − e^{−s}s − 2 e^{−s}}{(e^{−2 s} − 1) s^{2}}

Problem Notes on 110 – 119.

The methods employed in the solutions involve partial fractions, the shifting theorem, differentiation theorem, Heaviside theorem, and the convolution theorem.

Some of the answers below were obtained by hand manipulation of the maple answers. The complications of different answers can be resolved by modifying the following example:

readlib(laplace): J:=invlaplace:        # Read library code
J((s-1)/((s-1)^2-4),s,t):               # Compute inverse transform
subs(cosh(t)=exp(t)/2+exp(-t)/2,"):     # Substitute for hyperbolic cosine
p:=simplify(");                         # Simplify expression (maybe)
q:=(exp(3*t)+exp(-t))/2:                # Claimed answer
simplify(p-q);                          # Zero if p equals q

Answer to 110. \frac{3}{s^{2}} = 3L(t). Answer: 3 t

Answer to 111. \frac{s + 1}{s^{2} + 1} = \frac{s}{1 + s^{2}} + \frac{1}{1 + s^{2}} = L(cos(t)) + L(sin(t)). Answer: cos(t) + sin(t)

Answer to 112. \frac{1}{(s − 2)^{3}} = (\frac{1}{s^{3}} ) |_{s→(s−2)} = (L(t^{2}/2)) |_{s→(s−2)} , then apply the shifting theorem

on the last term. Answer: \frac{t^{2}e^{2 t}}{2}

Answer to 113. \frac{3s + 2}{s^{2} + 6s + 25} = \frac{3s + 2}{(s + 3)^{2} + 4^{2}} = \frac{3(s + 3) + (−7/4)(4)}{(s + 3)^{2} + 4^{2}} = 3\frac{s + 3}{(s + 3)^{2} + 4^{2}} +

(−7/4)\frac{4}{(s + 3)^{2} + 4^{2}} , then apply tables. Answer: 3 e^{−3 t} cos(4 t) − \frac{7 e^{−3 t} sin(4 t)}{4}

Answer to 114. \frac{s − 1}{(s − 1)^{2} − 4} is the shift of the transform of a hyperbolic function. Answer:

0.5e^{3t} + 0.5e^{−t}

Answer to 115. \frac{s + 1}{(s − 1)(s − 2)} = \frac{a}{s − 1} + \frac{b}{s − 2} = L(a exp(t) + b exp(2t)). To determine a and

b apply partial fraction theory. Answer: 3e^{2t} − 2e^{t}

Answer to 116. ln (1 + \frac{4}{s^{2}} ) uses a standard trick: let F (s) = ln(1 + 4s^{−2}), then F′(s) =

−8s^{−3}/(1 + 4s^{−2}) = −8/(s(s^{2} + 4)) = \frac{a}{s} + b\frac{s}{s^{2}+4} + c\frac{2}{s^{2}+4} for some a, b, c using the theory of partial fractions. Therefore, (d/ds)L(f(t)) = L(a + b cos(2t) + c sin(2t)) and by (d/ds)L(f(t)) =

L((−t)f(t)) it follows that (−t)f(t) = a + b cos(2t) + c sin(2t). Answer: −\frac{−2 + 2 cos(2 t)}{t}

Answer to 117. arctan(1/s) uses the method of (g). Answer: \frac{sin(t)}{t}

Answer to 118. \frac{1}{s} arctan(1/s) uses the convolution method. Let g(t) be the answer to (h) above:

L(g(t)) = arctan(1/s). Then (1/s) arctan(1/s) = L(1)L(g(t)) = L(∫_{0}^{t}(1)g(r)dr) so the answer is

∫_{0}^{t} g(r)dr. Answer: Si(t) = ∫_{0}^{t} sin(r)/r dr

Answer to 119. \frac{exp(−2s)}{s} is the transform of a Heaviside function. Answer: H (t − 2)

Problem Notes on 120 – 127.

The answers and maple notes appear below. The methods employed in the solutions involve partial fractions, the shifting theorem, differentiation theorem, Heaviside theorem, and the convolution theorem. Some of the problems require knowledge of extended transform tables.

Answer to 120. \frac{2e^{−s}}{s^{2} + 16} = L(\frac{H (t − 1) sin(4 t − 4)}{2})

Answer to 121. \frac{1}{s}(1 − exp(−2s))^{2} = L(1 − 2 H (t − 2) + H (t − 4)). The maple code is simple yet

tricky enough to cause a stall. It relies upon the expansion by hand of the term (1 − exp(−2s))^{2}.

      readlib(laplace): J:=invlaplace:        # Read library code
      J(1/s -2*exp(-2*s)/s + exp(-4*s)/s,s,t);

Answer to 122. \frac{1}{s} tanh(s) = L(f(t)) where f(t) is the 4-periodic extension of a square wave H(t) − 2H(t − 2) on [0, 4] to the whole real line. Some versions of maple return

invlaplace(\frac{tanh(s)}{s}, s, t)

which means no evaluation was done. It doesn’t mean it can’t be done, only that maple can’t figure it out. This is one of those problems that relies upon knowledge of transform theory and in particular intuition about the transform of a square wave gained by hand calculation or mathematical tables.

Answer to 123. \frac{exp(−2s)}{s^{2} − 1} = L(sinh(t − 2)H(t − 2))

Answer to 124. \frac{1 − exp(−4s)}{s^{2}} = L(t − (t − 4)H(t − 4)). The answer is checked by this maple code:

readlib(laplace):                       # Read library code
J:=invlaplace: L:=laplace:
readlib(Heaviside):                     # Read in library function
H:=x -> Heaviside(evalf(x)):            # Abbreviation
# Write as (1/s^2) - (1/s)*(exp(-4*s)/s) and use convolutions
p:= unapply(J(1/s,s,u),u);              # L(p)=1/s
q:= unapply(J(exp(-4*s)/s,s,u),u);      # L(q)=exp(-4*s)/s
J(1/s^2,s,t) - int(p(t-x)*q(x),x=0..t); # Convolution theorem
# Check our work
lprint(`F(s) = (1-exp(-4*s))/s^2`);     # The required transform F(s)
L(t,t,s)-L((t-4)*H(t-4),t,s);           # Maple's idea of F(s)

Answer to 125. \frac{s(1 − exp(−2πs))}{s^{2} + 1} = L(cos(t) − cos(t)H(t − 2π))

The solution is similar to the previous problem. The maple changes:

p:=unapply(J(s/(s^2+1),s,u),u);
q:=unapply(J(exp(-2*Pi*s)/s,s,u),u);
J(1/(s^2+1),s,t) - int(p(t-x)*q(x),x=0..t);

Answer to 126. \frac{1}{(s^{2} + 1)^{3}} + s^{−5} exp(−10s) = L(f(t)) where f(t) = −\frac{t^{2} sin(t)}{8} − \frac{3 t cos(t)}{8} +

\frac{3 sin(t)}{8} + \frac{H (t − 10) (t − 10)^{4}}{24}

Answer to 127. \frac{1}{s^{4}(s^{2} + 1)} = L(\frac{t^{3}}{6} − t + sin(t))

Problem Notes on 128 – 133.

The differential equation problems are worked by hand whenever possible. Maple is used to check answers and to do clumsy computations. Answer to 128. y′′ + 4y = 0, y(0) = 0, y′(0) = 2, y(t) = sin(2 t).

Answer to 129. y′′ − 4y = 0, y(0) = 1, y′(0) = 1, y(t) = (3/4)e^{2t} + (1/4)e^{−2t}.

Answer to 130. y′′ + 3y′ + 2y = 0, y(0) = 0, y′(0) = 1, y(t) = e^{−t} − e^{−2t}.

Answer to 131. y′′ + 3y′ + 2y = exp(−t), y(0) = y′(0) = 0, y(t) = te^{−t} − e^{−t} + e^{−2 t}.

Answer to 132. y′′ − 5y′ + 4y = exp(3t), y(0) = 1, y′(0) = −1, y(t) = \frac{11 e^{t}}{6} − \frac{e^{3 t}}{2} − \frac{e^{4 t}}{3} .

Answer to 133. y′′ + 2y′ + 2y = t + sin(t), y(0) = y′(0) = 0, y(t) = \frac{t}{2} − 1/2 + \frac{sin(t)}{5} − \frac{2 cos(t)}{5} + \frac{e^{−t} sin(t)}{5} + \frac{9 e^{−t} cos(t)}{10}.

A model example for dealing with differential equations appears below.

# Example for y'' - 4y=0, y(0)=1, y'(0)=1
readlib(laplace): J:=invlaplace: L:=laplace:    # Read library code
de:=diff(y(t),t,t)-4*y(t)=0:                    # Define differential equation
ic:=y(0)=1,D(y)(0)=1:                           # Define initial conditions
ans:=L(de,t,s);                                 # Transform the DE
ans:=subs({ic},ans);                            # Substitute the initial conditions
Ly:=solve(ans,L(y(t),t,s));                     # Solve for transform of y
J(Ly,s,t);                                      # Compute inverse transform

It is also possible to find only the answer and not the intermediate steps of the solution by the transform method. This method uses dsolve. It is recommended for checking answers obtained by other methods.

# Example for y'' - 4y=0, y(0)=1, y'(0)=1
de:=diff(y(t),t,t)-4*y(t)=0:                    # Define differential equation
ic:=y(0)=1,D(y)(0)=1:                           # Define initial conditions
dsolve({de,ic},y(t),laplace);                   # Solve by laplace method

Shown below is a method for checking other forms of the same answer.

# Example for y'' - 4y=0, y(0)=1, y'(0)=1
p:=2*(cosh(t))^2 - 1 + sinh(t)*cosh(t):         # First answer
q:=(3*exp(2*t)+exp(-2*t))/4:                    # Claimed answer
simplify(p-q);                                  # Zero if p equals q (sometimes fails)
subs(cosh(t)=exp(t)/2+exp(-t)/2,p):             # Substitute for cosh(t)
r:=subs(sinh(t)=exp(t)/2-exp(-t)/2,"):          #  and sinh(t)
simplify(r-q);                                  # Zero if r equals q (sometimes fails)
plot(r-q,t=0..1);                               # Test theory with a plot

In some versions of maple the above procedure never resolves the question of equality of the two answers p and q. Even the plot is strange (it should be the zero plot). However, to be fair, it is possible to look at the maple output and decide it is zero, i.e., p = q.

Problem Notes on 134 – 139.

The transformed equation is the function of s which equals L(y(t)). The various methods used below involve differential equations, Heaviside relations and the standard theorems from transform theory.

Answer to 134. To solve Y′ + (1/(s − 1))Y = 0, Y (0) = 2 compute the integrating factor exp(∫ 1/(s − 1)ds) = s − 1, multiply by it, then factorize to obtain ((s − 1)Y)′ = 0, L(y) = Y = −2/(s − 1), y(t) = −2 exp(t). To check the answer use these maple steps:

      readlib(laplace): J:=invlaplace: L:=laplace:    # Read library code
      L(-2*exp(t),t,s);                               # Transform claimed answer
      Y:=unapply(",s);                                # Make transform a function
      diff(Y(s),s)+(1/(s-1))*Y(s); Y(0);              # Check DE and IC

Answer to 135. To solve the equation ty′′ − ty′ − y = t, y(0) = 0, y′(0) = 1 apply the formula L((−t)f(t)) = (d/ds)L(f(t)) to make the following simplifications of the transformed equation:

(−d/ds)L(y′′) + (d/ds)L(y′) − L(y) = 1/s^{2}

(d/ds)[−(s(sL(y) − y(0)) − y′(0)) + sL(y) − y(0)] − L(y) = 1/s^{2}

Y = L(y), Y′ + (−2/(1 − s))Y = 1/(s^{3}(1 − s))

The differential equation can be solved by hand. The answer is

Y (s) = L(y) = \frac{(1/2)(2s − 1)}{(s^{2})(s − 1)^{2}} + \frac{a}{(s − 1)^{2}}

where a is an arbitrary constant. Partial fraction theory can be applied to find y. The result is an answer for y containing the parameter a. Using the initial conditions y(0) = 0 and y′(0) = 1 results in a = 1, therefore the solution is given by

y(t) = (3/2)t exp(t) − (1/2)t

The basic ideas for using maple appear in the example below.

# A similar example.
# ty'' -2ty' - 2y = 8t, y(0)=0, y'(0)=2
# Y(s)=L(y), Y' + (2/(s-2))Y=(-8)/(s^3(s-2)) [done by hand]
readlib(laplace): J:=invlaplace: L:=laplace:        # Read library code
Y:='Y': y:='y': a:='a':                             # Reset variables
de:=diff(Y(s),s)+(2/(s-2))*Y(s)=(-8)/((s^3)*(s-2)): # Define DE
dsolve(de,Y(s)):                                    # Solve DE
subs(_C1=a,"):                                      # Substitute for arb const
Y:=unapply(rhs("),s);                               # Make answer a function
# Y(s) = (8*s-8)/((s^2)*(s-2)^2) + a/(s-2)^2
J(Y(s),s,t);                                        # Find inverse transform
y:=unapply(",t);                                    # Make answer a function
# y := t -> 2*t*exp(t) + t*exp(t)* a - 2*t:
simplify(t*diff(y(t),t,t)-2*t*diff(y(t),t)-2*y(t) - 8*t); # Check the DE
y(0); D(y)(0);                                      # Print the IC
# Choose a=2 to make the IC work: y(0)=0, y'(0)=2
# y(t) = 4*t*exp(t) - 2*t

The partial fraction problem can be solved in maple as follows

Y:= s -> (1/2)*(2*s-1)/((s^2)*(s-1)^2) + a/(s-1)^2: # Define function of s
convert(Y(s),parfrac,s);                            # Convert to partial fractions

Answer to 136. The equation to solve is y′′ + 6y′ + 9y = 10 sin(10t), y(0) = 0, y′(0) = 0. The transform Y (s) = L(y(t)) is given by Y = 100/((s^{2} + 100)(s + 3)^{2}). The answer is

y(t) = \frac{100 te^{−3 t}}{109} + \frac{600 e^{−3 t}}{11881} − \frac{910 sin(10 t)}{11881} − \frac{600 cos(10 t)}{11881}

Answer to 137. The transform of y′′′ + y = 0, y(0) = y′(0) = 0, y′′(0) = 1 is Y (s) = L(y(t)) = 1/(s^{3} + 1) = 1/((s + 1)((s − 1/2)^{2} + 3/4)).

y(t) = \frac{e^{−t}}{3} + \frac{e^{\frac{t}{2}} sin(\frac{\sqrt{3}t}{2})\sqrt{3}}{3} − \frac{e^{\frac{t}{2}} cos(\frac{\sqrt{3}t}{2})}{3}

Answer to 138. The transform of y^{(4)} − y = 0, y(0) = 0, y′(0) = y′′(0) = 0, y′′′(0) = −1 is Y (s) = L(y(t)) = −1/(s^{4} − 1) = −1/((s − 1)(s + 1)(s^{2} + 1). Partial fractions can be applied to find y(t) = (−1/4) exp(t) + (1/4) exp(−t) + (1/2) sin(t).

Answer to 139. The transform of y^{(4)} − y = cos(t), y(0) = y′′(0) = 1, y′(0) = y′′′(0) = 0 is Y (s) = L(y(t)) = (s^{3} + 3 + s/(s^{2} + 1))/(s^{4} − 1). The answer is

y(t) = (5/8) exp(t) + (5/8) exp(−t) − (1/4) cos(t) − (1/4)t sin(t)

Problem Notes on 140 – 145.

The solutions to the problems use standard transform theorems and formulas, the formula for the transform of a periodic function, and the convolution theorem.

Answer to 140. The right side of the given formula is a direct application of the theory of the transform of a periodic function: L(f(t)) = ∫_{0}^{T} e^{−st}f(t)dt/(1 − exp(−T s). The function f(t) is a periodic function of period T = π built from the curve sin(t) on 0 ≤ t ≤ T, i.e., f(t) = | sin(t)|. The left side can be written as L(−ty(t)). The answer is y(t) = −| sin(t)|/t.

Answer to 141. The differential equation (d/ds)^{2}L(y) = F (s) is equivalent to L((−t)^{2}y(t)) = F (s). Therefore, y(t) = (1/t^{2})f(t) where L(f(t)) = F (s). Applied to

F (s) = \frac{2 e^{−2 s−2}}{s + 1} + \frac{e^{−2 s−2}}{(s + 1)^{2}}

gives y(t) = exp(−t)H(t − 2)/t.

Answer to 142. The differential equation y′′ + y = H(t − 10), y(0) = y′(0) = 0 transforms to Y (s) = L(y) = (1/(s^{2} + 1)) · L(H(t − 10)). The answer can be given by the convolution theorem because 1/(s^{2} + 1) = L(sin(t)). The answer is y(t) = (1 − cos(t − 10))H(t − 10).

Answer to 143. The equation y′′ + 10y′ + 50y = 25t(1 − H(t − 4)) + 50H(t − 4), y(0) = y′(0) = 0 is solved by the method of convolutions. The answer is y(t) = ∫_{0}^{t} k(t − r)f(r)dr where L(k(t)) = 1/(s^{2} + 10s + 50) and f(t) = 25t + (50 − 25t)H(t − 4). After integration

y(t) = (\frac{11}{10} + e^{20−5 t} sin(5 t − 20) + \frac{9 e^{20−5 t} cos(5 t − 20)}{10} − \frac{t}{2}) H(t−4)+ \frac{t}{2} −1/10+ \frac{e^{−5 t} cos(5 t)}{10}.

Answer to 144. The idea for solving the equation

y(t) = exp(−t) − ∫_{0}^{t} sin(t − r)y(r)dr

is to apply L to both sides. The right side contains a term that is transformed by the convolution theorem. After some manipulations

L(y) = \frac{1}{s + 1} − \frac{1}{s^{2} + 1}L(y)

y(t) = \frac{2 e^{−t}}{3} + \frac{cos(\sqrt{2}t)}{3} − \frac{sin(\sqrt{2}t)\sqrt{2}}{6}

Answer to 145. The equation

y(t) = −t + ∫_{0}^{t} \frac{(t − r)^{3}}{6} y(r)dr

is solved by the same method as 7(e) above:

L(y) = \frac{−1}{s^{2}} + \frac{1}{s^{4}} L(y)

y(t) = \frac{e^{−t}}{4} − \frac{e^{t}}{4} − \frac{sin(t)}{2}

Problem Notes on 146 – 151.

To solve for y below means a manipulation is required to arrange the given equation in the form L(y) = L(f(t)). Lerch’s theorem is applied to conclude y = f(t). In each case the theory of partial fractions is required. Answer to 146. The theory of partial fractions says that

L(y) = \frac{s + 1}{(s^{2} + 1)(s + 3)^{2}} = \frac{as + b}{s^{2} + 1} + \frac{c}{s + 3} + \frac{d}{(s + 3)^{3}}

Various hand manipulations can be used to determine a, b, c, d. Maple can be used to check hand calculations or to do the initial decomposition. Answer: a = 1/50, b = 7/50, c = −1/50, d = −1/5,

y(t) = \frac{cos(t)}{50} + \frac{7 sin(t)}{50} + \frac{−e^{−3 t}}{50} + \frac{−te^{−3 t}}{5}

Using maple to decompose a fraction into partial fractions followed by solving for y:

      p:=(s+1)/((s^2+1)*(s+3)^2) :                  # Fraction to decompose
      pf:=convert(p,parfrac,s);                     # Decomposition with numbers
      readlib(laplace): L:=laplace: J:=invlaplace:  # Read library code
      J(pf,s,t);                                    # Solve L(y)=pf for y

Using maple to check a hand computation:

      p:=(s+1)/((s^2+1)*(s+3)^2) :                    # Fraction to decompose
      q:=(a*s+b)/(s^2+1)+c/(s+3)+d/(s+3)^2:           # Decomposition with symbols
      subs({a=1/50,b=7/50,c=-1/50,d=-1/5},expand(q)); # Put answers into fractions
      simplify(");                                    # Simplify fractions
      simplify("-p);                                  # Expect zero

Answer to 147. The answers:

L(y) = \frac{s}{(s^{2} + 1)(s + 1)^{2}} = \frac{1/2}{s^{2} + 2} + \frac{−1/2}{(s + 1)^{2}}

y(t) = \frac{sin(t)}{2} − \frac{te^{−t}}{2}

Answer to 148. The answers:

L(y) = \frac{2}{(s + 1)(s − 1)(s + 2)} = \frac{−1}{s + 1} + \frac{1/3}{s − 1} + \frac{2/3}{s + 2}

y(t) = −e^{−t} + \frac{1}{3}e^{t} + \frac{2}{3}e^{−2t}

Answer to 149. The answers:

L(y) = \frac{exp(−4s)}{(s − 1)(s + 1)(s − 2)} = e^{−4 s} (\frac{1}{6 s + 6} − \frac{1}{2 s − 2} + \frac{1}{3 s − 6})

y(t) = (\frac{e^{4−t}}{6} − \frac{e^{t−4}}{2} + \frac{e^{2 t−8}}{3} ) H(t − 4)

Answer to 150. Transform L(y) + L(t sin(t)) = F (s) where F (s) = \frac{s+1}{(s−1)(s+2)^{2}(s^{2}+1)} into L(y + t sin(t)) = L(f) where L(f) = F (s). The fraction decomposition of F is

F (s) = \frac{1/9}{s − 1} + \frac{1/15}{(s + 2)^{2}} + \frac{2/225}{s + 2} − \frac{4/25 + (3/25)s}{s^{2} + 1}

The answer is obtained from y(t) + t sin(t) = f(t):

y(t) = \frac{e^{t}}{9} + \frac{te^{−2 t}}{15} + \frac{2 e^{−2 t}}{225} − \frac{4 sin(t)}{25} − \frac{3 cos(t)}{25} − t sin(t)

Answer to 151. The answers:

L(y) = (\frac{6}{s^{3}} − \frac{2}{s + 1}) \frac{s − 1}{s} = \frac{6}{s^{3}} − \frac{6}{s^{4}} + \frac{2}{s} − \frac{4}{s + 1}

y(t) = 3 t^{2} − t^{3} + 2 − 4 e^{−t}

Problem Notes on 152 – 155.

To solve for y below involves arrangement of the given equation into the form L(y) = L(u(t)). Lerch’s theorem is applied to conclude y = u(t). The differentiation theorem, the theory of partial fractions and the convolution theorem are the main tools.

Answer to 152. The problem L(y) = \frac{s}{(s−1)(s+2)}L(tH(t − 1)) can be solved by the convolution theorem as y(t) = ∫_{0}^{t} f(t − r)g(r)dr provided L(f) = \frac{s}{(s−1)(s+2)} and L(g) = L(tH(t − 1)). The answer:

y(t) = (−1/2 + \frac{2}{3}e^{t−1} − \frac{1}{6}e^{2−2t}) H(t − 1)

Answer to 153. The problem L(y′′(t))+sy(0)+y′(0) = \frac{s}{s+1} is equivalent to L(y) = 1/(s(s+1)). The answer: y(t) = 1 − exp(−t).

Answer to 154. The differential equation \frac{d}{ds}L(−y) + \frac{d^{2}}{ds^{2}} L(y) = F (s) has constant coefficients and therefore it can be rearranged in the form L((−t + t^{2})y(t)) = F (s). The answer is y(t) = f(t)/(t^{2} + t) where L(f(t)) = F (s) = \frac{s−1}{(s+1)(s^{2}+2s+2)}. The answer:

y(t) = \frac{e^{−t} (sin(t) − 2 + 2 cos(t))}{t (t + 1)}

Answer to 155. The problem L(y) = L(t exp(−t))\frac{1}{s} arctan(s) is an application of the convolution solution y(t) = ∫_{0}^{t} f(t−r)g(r)dr where L(f) = L(t exp(−t)) and L(g) = arctan(s)/s. The answer:

y(t) = ∫_{0}^{t}−Si(t − r)re^{−r}dr

where Si(x) = ∫_{0}^{x} sin(t)/tdt.

Problem Notes on 156 – 159.

The differentiation theorem, the theory of partial fractions and the convolution theorem are the main tools from transform theory. Linear algebra (or Cramer’s Rule) is required to solve the linear systems of algebraic equations.

Answer to 156. The integral equation L(∫_{0}^{t} sin(t − r)y(r)dr) = \frac{s}{(s^{2}+1)^{2}} can be simplified by the convolution theorem to give L(y) = \frac{s}{s^{2}+1}. The answer is obtained from tables.

Answer to 157. The system of equations L(x) + sL(y) = 1, (s − 1)L(x) − sL(y) = 0 can be solved by Cramer’s Rule or Gaussian elimination to give the transform answers L(x) = 1/s, L(y) = (s − 1)/s^{2}. Inversion of these formulas gives the answers x(t) = 1, y(t) = 1 − t.

Answer to 158. The equation L(∫_{0}^{t} y(r)dr) + 20L(y) = F (s) is equivalent by the convolution theorem to L(1)L(y) + 20L(y) = F (s) and therefore

L(y) = (\frac{1}{s^{2}} − \frac{exp(−s)}{s^{2}} − \frac{exp(−s)}{s})/(20 + 1/s)

The solution by inversion is

y(t) = 1 − e^{−t/20} − H(t − 1) + \frac{19}{20}e^{(1−t)/20})H(t − 1)

Answer to 159. The system of equations L(x) + L(y) = \frac{s}{(s−1)(s+1)}, L(x) − L(y) = \frac{s}{(s−2)(s+2)} can be solved by Gaussian elimination to give the transform answers

L(x) = \frac{s/2}{(s − 1)(s + 1)} + \frac{s/2}{s − 2)(s + 2)}

L(y) = \frac{s/2}{(s − 1)(s + 1)} − \frac{s/2}{s − 2)(s + 2)}

and therefore the answers are

x(t) = \frac{cosh(t)}{2} + cosh^{2}(t) − 1/2

y(t) = \frac{cosh(t)}{2} − cosh^{2}(t) + 1/2

Problem Notes on 160.

The mathematical analysis is carried out by hand. The use of maple is possible, however it is somewhat more complicated.

The purpose of this conversion is to put the differential equations into standard form. Almost every symbolic and numeric solver requires the standard form before any solution can be done. Maple is no exception: the standard form is required for input into dsolve.

Problem Notes on 161.

The solution of the system of equations can be carried out by hand. Maple is recommended because of the complexity of the solution. Below is a maple program to compute the solution x(t), y(t).

L1:=1; L2:=1; R1:=1; R2:=1/4; E:=120; M:=1/2;   # Constants
K:=L1*L2-M*M;                                   # Definition of K,
a:= -R1*L2/K; b:=M*R2/K; f:=E*L2/K;             #  a, b, f
c:=M*R1/K; d:= -R2*L1/K; g:= -M*E/K;            #  c, d, g
de1:=diff(x(t),t)=a*x(t)+b*y(t)+f:              # First DE
de2:=diff(y(t),t)=c*x(t)+d*y(t)+g:              # Second DE
ic:=x(0)=0,y(0)=0:                              # Initial conditions
dsolve({de1,de2,ic}, {x(t),y(t)}, laplace);     # Solve for x, y

Problem Notes on 162.

The primary current is x(t) = I_{1}(t) and the secondary current is y(t) = I_{2}(t). The steady state value of a current is its limit at t = ∞. The question can be rephrased as follows: find t_{1} (about 9) and t_{2} (about 22) such that

0.99(120) = x(t_{1})
0.01(−39.23) = y(t_{2})

These equations can be solved experimentally with a graph or by using sophisticated root–solving algorithms. The fsolve routine in maple can do both problems. From the graph of y it is seen that t_{2} > 1.115. The number 1.115 is an approximate t-value for the maximum excursion of y (min y(t) is about −39.23). A more precise value can be obtained by solving y′(t) = 0 for t: tanh(\sqrt{13}t/6) = 80\sqrt{13}/400.

Problem Notes on 163.

To show that V (t) satisfies the differential equation divide the given differential equation by C. To show that V (t) satisfies the initial conditions V (0) = V′(0) = 0 apply the given initial conditions.

Problem Notes on 164.

The objective is to find k in the equation L(k(t)) = \frac{1}{s^{2} + 2s + 2} by inverse transform theory. The

answer is k(t) = e^{−t} sin(t).

Problem Notes on 165.

The convolution integral

V (t) = \frac{1}{C} ∫_{0}^{t} k(t − x)E(x)dx

has transform L(V) = L(k)L(E/C) or [Ls^{2} + Rs + 1/C]L(V) = L(E)/C. Obtain the latter by transforming the differential equations and initial conditions and then reverse the steps.

Problem Notes on 166.

The problem reduces to integration of

V (t) = 20 ∫_{0}^{t} exp(x − t) sin(t − x) exp(−x) sin(5x)dx

The answer is

V (t) = \frac{10 sin(t)}{3 e^{t}} − \frac{10 cos^{2}(t) sin(t) (4 cos^{2}(t) − 3)}{3 e^{t}}

or

V (t) = \frac{10}{3} sin^{3}(t)e^{−t}(1 + 4 cos^{2}(t)).

Problem Notes on 167.

This problem is one that should not be attempted by hand. The technical difficulties associated with hand calculation and graphing are substantial.

A zero of y is a point x where y(x) = 0. A critical point of y is a point x where y′(x) = 0 (generally a maximum or minimum of y). The sample code below implements the Heun Method on a similar problem. To implement RK4 replace heun by rk4 in the code.

# Solve  y' = y*y + g(x), y(0)=0, where g(x) is a given function.
# Use the Heun Method step size h on [a,b].
#
read `/u/cl/maple/desolvers`:                   # Define Heun and RK4 methods
g:= x -> -cos(Pi*x)*Pi+2*cos(2*Pi*x)*Pi:        # Special function g(x)
F:= (t,y) -> y^2 + g(t):                        # Right side of DE
a:=0.0: b:=1.0: m:=100: h:=evalf((b-a)/m):      # Solve on [a,b] step h
x0:=a: y0:=0.0:                                 # Initial conditions
L:=[x0,y0]:                                     # Make point list
for i from 1 to m do
  y0:=evalf(heun(x0,y0,h,F)): x0:=evalf(x0+h): L:=L,[x0,y0]:
od:
plot([ L ]);                                    # plot the numeric solution
print([ L ]);                                   # Find zeros and critical points

The source code for Heun and RK4 appears below. A machine readable copy is in the file /u/cl/maple/desolvers. Detailed study of Euler’s method and Heun’s method is worthwhile because intuition can be gained. A deep understanding of the Runge-Kutta method is possible with additional references or a specialized course in numerical analysis.

# Heun's Method
# usage: f:=proc(x,y) a maple procedure for F(x,y)
#        h=step, (x0,y0)=initial point.
heun:=proc(x0,y0,h,f)
local s,s0:
 s0 := f(x0,y0):

 s  := y0+h*(s0+f((x0+h),(y0+h*s0)))/2:
 RETURN(s):
end:
#
# Runge-Kutta Method order 4.
# usage: f:=proc(x,y) a maple procedure for F(x,y)
#        h=step, (x0,y0)=initial point.
rk4:=proc(x0,y0,h,f)
local k1,k2,k3,k4,s:
 k1:=h*f(x0,y0):
 k2:=h*f((x0+h/2),(y0+k1/2)):
 k3:=h*f((x0+h/2),(y0+k2/2)):
 k4:=h*f((x0+h),(y0+k3)):
 s:=y0+(k1+2*k2+2*k3+k4)/6:
 RETURN(s):
end:

Tables printed in maple can have bad column alignment. The function printf is used to produce formatted output with column alignment. Maple sources that use print can be modified to use printf, which C and Pascal programmers will find familiar. In some cases it is less work to save the table in a file and edit it in emacs to obtain the column alignment. Using print a list of plot points can be printed as a table as follows:

LL:=[ L ]:                      # Make an array from list L
N:=nops(LL):                    # Get size of list
for i from 1 to N do
  print(LL[i][1],LL[i][2]);     # LL[i] is a list of the form [x,y]
od:                             #   and LL[i][1] is x, LL[i][2] is y

Problem Notes on 168.

This problem is one that should not be attempted by hand. The technical difficulties associated with hand calculation and graphing are substantial.

The linear interpolant g of the data points (0, 3), (1/3, −5), (2/3, −6), (1, 10) is the function defined by the equation

g(t) = \{ \begin{matrix}−24t + 3 & 0 ≤ t < 1/3 \\ −3t − 4 & 1/3 ≤ t < 2/3 \\ 48t − 38 & 2/3 ≤ t ≤ 1 \end{matrix}

The function g goes through the data points and is linear on each subinterval. Heaviside’s unit step function H is used to express g in a different form. The method depends upon the representation

H(t − a) − H(t − b) = \{ \begin{matrix}1 & a ≤ t < b \\ 0 & t < a, t ≥ b \end{matrix}

If

g(t) = \{ \begin{matrix}f_{1}(t) & t_{0} ≤ t < t_{1} \\ f_{2}(t) & t_{1} ≤ t < t_{2} \\ f_{3}(t) & t_{2} ≤ t ≤ t_{3} \end{matrix}

then on t_{0} ≤ t ≤ t_{3}

g(t) = f_{1}(t)[H(t − t_{0}) − H(t − t_{1})] + f_{2}(t)[H(t − t_{1}) − H(t − t_{2})] + f_{3}(t)[H(t − t_{2}) − H(t − t_{3})] + f_{3}(t)H(t − t_{3})

This form is equivalent on [t_{0}, t_{3}] to the simpler form

g(t) = f_{1}(t)H(t − t_{0}) + [f_{2}(t) − f_{1}(t)]H(t − t_{1}) + [f_{3}(t) − f_{2}(t)]H(t − t_{2})

which becomes for the given data the function

g(x) = (−24t + 3)H(t) + (21t − 7)H(t − 1/3) + (51t − 34)H(t − 2/3)

The method is extensible to more than three intervals. It is computationally efficient to apply the method on five intervals or less. For a general method of handling data sets of arbitrary size see the solutions below.

The numerical solution and graphing can be done with a hand calculator, although this is certain to be tedious. The calculation by hand produces a table of values for plotting, from which the maximum and mimimum values of y can be determined.

An implementation appears below for maple on a similar problem using the same kind of function g. It is recommended to code the example, get it to work, then modify the source to do the problem. To locate the maximum and minimum from the maple solution it is necessary to print the list L and search for the largest and smallest y-values. The method was illustrated in aprevious problem.

a:=0.0: b:=1.0:                                 # Define [a,b]
readlib(Heaviside):                             # Read in library function Heaviside
H:=x -> Heaviside(evalf(x)):                    # Abbreviation
g:=t-> (-24*t + 3)*H(t) + (21*t-7)*H(t-1/3) +   # Define experimental function g
       (53*t-34)*H(t-2/3):                      #  in terms of Heaviside function
F:= (t,y) -> y^2 + g(t):                        # Define right side of DE
read `/u/cl/maple/desolvers`:                   # Get RK4 source code from library
m:=50: h:=evalf((b-a)/m):                       # Set number of steps and step size h
x0:=a: y0:=0.0:                                 # Initial data for DE
L:=[x0,y0]:                                     # Initial plotlist
for i from 1 to m do
  y0:=evalf(rk4(x0,y0,h,F)): x0:=evalf(x0+h):   # Apply RK4, compute (x0,y0)
  L:=L,[x0,y0]:                                 # Add (x0,y0) to plot list
od:
plot([ L ]);                                    # plot numeric solution (brackets required)

Problem Notes on 169.

This problem is one that should not be attempted by hand. The technical difficulties associated with hand calculation and graphing are substantial.

An implementation appears below for maple on a similar problem using the same kind of function g. The run time is about 15 seconds. It is recommended to code the example, get it to work, then modify the source to do the problem. To locate the maximum and minimum it is necessary to print the list L and search for the largest and smallest y-values. Computing wizards might modify the loop below to find the maximum and minimum in a more elegant fashion.

n:=12: a:=0.0: b:=1.0:                          # Define [a,b], n
X:=array(0..n,[seq(i*(b-a)/n,i=0..n)]):         # x-values

Y:=array(0..n,[                                 # y-values
  3.141592654,   2.408784785,  0.4822403291,    # Arrays are entered
 -0.7319568225, -3.224726277, -5.360718345,     #  in maple with this
 -6.807829539,  -7.298315887, -6.480432161,     #  special syntax: dimension
 -4.051443799,  -0.1464077496, 4.361205273,     #  first and then a list
  6.393106333                                   #  of values in brackets [].
]):
readlib(Heaviside):                             # Read in library function Heaviside
H:=x -> Heaviside(evalf(x)):                    # Abbreviation
F:= (t,y) -> y^2 + g(t):                        # Right side of DE
whole:= x -> convert(x,fraction):               # Convert decimal to integer
K1:= t -> whole(sum(H(t-X[kk]),kk=1..n-1)):     # K1 is the staircase function
P1:= (t,k) -> evalf((Y[k+1]-Y[k])*(t-X[k])/     # P1 is a straight line
                    (X[k+1]-X[k])+Y[k]):        #  (a linear function)
g:= t ->  P1(t,K1(t)):                          # Linear interpolant
read `/u/cl/maple/desolvers`:                   # Get RK4 source code
m:=50: h:=evalf((b-a)/m):                       # Set step size
x0:=a: y0:=0.0:                                 # Initial data for DE
L:=[x0,y0]:                                     # Initial plotlist
for i from 1 to m do
  y0:=evalf(rk4(x0,y0,h,F)): x0:=evalf(x0+h):   # Apply RK4, compute (x0,y0)
  L:=L,[x0,y0]:                                 # Add (x0,y0) to plot list
od:
plot([ L ]);                                    # plot numeric solution (brackets required)

Problem Notes on 170.

The conversion from scalar equation to system is always required for numerical procedures. There is a general procedure that applies to differential equations of order n: Let x_{i} = y^{(i)}(t) for 0 ≤ i ≤ n − 1, then y^{(n)}(t) = f(t, y(t), y′(t), …, y^{(n−1)}(t)) is equivalent to the system

x′_{0}  = x_{1}
x′_{1}  = x_{2}
⋮  ⋮
x′_{n−2} = x_{n−1}
x′_{n−1} = f(t, x_{0}, x_{1}, …, x_{n−1})

For order two, x_{0} = y(t) is the position and x_{1} = y′(t) is the velocity, therefore the procedure is sometimes called the position–velocity substitution.

Problem Notes on 171.

By problem 170, this system is equivalent to a second order scalar equation with constant coefficients. The answers x and y are the position u(t) and velocity u′(t) for the second order scalar equation u′′ + 2u′ + 3u = exp(t), u(0) = 1, u′(0) = 2.

To solve the system

x′ = y, y′ = −2x − 3y + exp(t), x(0) = 1, y(0) = 2

by the Runge-Kutta method for x(1), y(1) using step size 0.2 it is convenient to use a maple library for numeric solutions. Illustrated below is a slightly different equation.

# Solve x'=y, y'=2x+3y+sin(t), x(0)=1/5, y(0)=3/5 by RK4
read `/u/cl/maple/desolvers`:           # Read library code "rk4phase"

F:=(t,w) -> w[2]:                       # Right side of first DE
G:=(t,w) -> 2*w[1]+3*w[2] +sin(t):      # Right side of second DE
z:=array(1..2):                         # Array for answers
z[1]:=1/5: z[2]:=3/5:                   # Initial data
t0:=0:                                  # Initial point
n:=5:                                   # Use 5 steps
a:=0: b:=1:                             # Solve from t=a to t=b
h:=evalf((b-a)/n):                      # Use n steps of size h
# Compute numeric solutions x(t), y(t)
 L:=[t0,z[1],z[2]]:                     # Build list of answers
 for i from 1 to n do
   rk4phase(t0,z,h,F,G):                # RK4 step, answers in z[1], z[2]
   t0:=evalf(t0+h):                     # New t-value
   L:=L,[t0,z[1],z[2]]:                 # Add point to list L
 od:
 LL:=[ L ]:                             # List of answers
N:=nops(LL):                            # Number of answers to print
for i from 1 to N do                    # Print a table with N rows
  print(LL[i][1], LL[i][2], LL[i][3]);  # Print row i of table
od:
#                     .0000000000, .2000000000, 0.600000000
#                     .2000000000, .3777841124, 1.264325549
#                     .4000000000, .7534259653, 2.666967804
#                     .6000000000, 1.541168434, 5.565152275
#                     .8000000000, 3.175621765, 11.50018832
#                     1.000000000, 6.540082231, 23.61086390

The answers for the present problem are

\begin{matrix}0 & 1 & 2 \\ .2000000000 & 1.281391938 & .9134868419 \\ .4000000000 & 1.396260246 & .2995322557 \\ .6000000000 & 1.421050902 & −.01017506904 \\ .8000000000 & 1.404917732 & −.1243820161 \\ 1.000000000 & 1.379499054 & −.1122432701 \end{matrix}

Problem Notes on 172. The system of circuit equations can be written as

x′ = −20x + 5y + 120, y′ = 5x − 20y, x(0) = 0, y(0) = 1

To apply Runge-Kutta and graph x and y for 0 ≤ t ≤ 2 apply a modification of the following example:

# Solve x'= -x + y +100, y'=x-y, x(0)=0, y(0)=1 by RK4
# Plot x, y on [0,1] using one set of axes.
#
read `/u/cl/maple/desolvers`:           # Read library code "rk4phase"
F:=(t,w) -> -w[1]+w[2]+100:             # Right side of first DE
G:=(t,w) -> w[1]-w[2]:                  # Right side of second DE
z:=array(1..2):                         # Array for answers
z[1]:=0: z[2]:=1:                       # Initial data
t0:=0:                                  # Initial point
n:=50:                                  # Use 50 steps
a:=0: b:=1:                             # Solve from t=a to t=b
h:=evalf((b-a)/n):                      # Use n steps of size h

# Compute numeric solutions x(t), y(t)
 L1:=[t0,z[1]]:                         # Build two lists of answers
 L2:=[t0,z[2]]:                         #  one for x, one for y
 for i from 1 to n do
   rk4phase(t0,z,h,F,G):                # RK4 step, answers in z[1], z[2]
   t0:=evalf(t0+h):                     # New t-value
   L1:=L1,[t0,z[1]]:                    # Add point to list L1
   L2:=L2,[t0,z[2]]:                    # Add point to list L2
 od:
plot({[ L1 ], [ L2 ]});                 # Plot solutions in one figure
# Exact solutions to this example:
#           x := t -> 50*t - 49/2 + (51/2)*exp(- 2*t)
#           y := t -> 50*t + 51/2 - (51/2)*exp(- 2*t)

Problem Notes on 173. To solve by Runge-Kutta and graph on [0, π] the solution to the oscillator problem

u′′(t) + 4u(t) = g(t), u(0) = u′(0) = 0

it is necessary to convert to the system

x′ = y, y′ = g(t) − 4x, x(0) = y(0) = 0

The problem is similar to the previous one, with the following modification to define the external force

g(t) = sin(t)H(t) + (cos(t) + 1 − sin(t))H(t − π/2)

readlib(Heaviside):                     # Read in library function Heaviside
H:=x -> Heaviside(evalf(x)):            # Abbreviation
g:= t -> sin(t)*H(t) +                  # Define external
         (cos(t)+1-sin(t))*H(t-Pi/2):   #  force g(t)
F:=(t,w) -> w[2]:                       # Right side of first DE
G:=(t,w) -> -4*w[1]+g(t):               # Right side of second DE
z:=array(1..2):                         # Array for answers
z[1]:=0: z[2]:=0:                       # Initial data
a:=0: b:=Pi:                            # Define [a,b] for plot

Problem Notes on 174.

The Airy Differential Equation u′′ + tu = 0 with initial conditions u(0) = 0, u′(0) = 1 converts into a system of first order differential equations by the transformation x = u(t), y = u′(t). The system is x′ = y, y′ = −tx, x(0) = 0, y(0) = 1.

Problem Notes on 175.

The system

u′′′ − t^{2}u′′ + (u′)^{2}u = sin(t), u(0) = u′(0) = 0, u′′(0) = 1,

converts into a system of three first order differential equations by the transformation x = u(t), y = u′(t), z = u′′(t). The initial conditions are x(0) = y(0) = 0, z(0) = 1. The differential equations are x′ = y, y′ = z, z′ = t^{2}z − xy^{2} + sin(t).

Problem Notes on 176.

The source for RK4 for systems is rk4phase in file /u/cl/maple/desolvers. The application to the Airy differential equation uses the system x′ = y, y′ = −tx, x(0) = 0, y(0) = 1. To produce the three-column table it is necessary to print only certain rows of a much larger table. Here is an example for the equation x′ = y, y′ = x + y, x(0) = 2, y(0) = −2 on [0, 3]:

read `/u/cl/maple/desolvers`:           # Read library code "rk4phase"
F:=(t,w) -> w[2]:                       # Right side of first DE
G:=(t,w) -> w[1]+w[2]:                  # Right side of second DE
z:=array(1..2):                         # Array for answers
z[1]:=2: z[2]:= -2:                     # Initial data
t0:=0:                                  # Initial point
n:=60:                                  # Use 60 steps
a:=0: b:=3:                             # Solve from t=a to t=b
h:=evalf((b-a)/n):                      # Use n steps of size h
# Compute numeric solution t, x(t), y(t)
 L:=[t0,z[1],z[2]]:                     # Build list of answers
 for i from 1 to n do
   rk4phase(t0,z,h,F,G):                # RK4 step, answers in z[1], z[2]
   t0:=evalf(t0+h):                     # New t-value
   L:=L,[t0,z[1],z[2]]:                 # Add point to list L
 od:
 LL:=[ L ]:                             # List of answers
N:=nops(LL):                            # Number of objects in list LL
W:=[0,0.4,0.8,1.2,1.6,2.0,2.4,2.8,3.0]: # t-values for printing
M:=nops(W):                             # Number of answers to print
for i from 1 to N do                    # Print a table with certain rows
 for j from 1 to M do                   # Search small table for values,
  if evalf(LL[i][1])=evalf(W[j]) then   #  and if it matches, then
   print(LL[i][1], LL[i][2], LL[i][3]); #  print row i of table.
  fi:
 od:
od: