Applications related to ordinary and partial differential equations

Martha L. Abell , James P. Braselton , in Mathematica by Example (Sixth Edition), 2022

6.2.4 Variation of parameters

Let S = { y 1 , y 2 } be a fundamental set of solutions for Eq. (6.9). To solve the nonhomogeneous equation (6.8), we need to find a particular solution, y p of Eq. (6.8). We search for a particular solution of the form

A particular solution, y p , is a solution that does not contain any arbitrary constants.

Key concept: A general solution of the nonhomogeneous linear equation is y = y h + y p , where y h is a general solution of the corresponding homogeneous equation, and y p is a particular solution of the nonhomogeneous equation.

(6.14) y p = u 1 ( t ) y 1 ( t ) + u 2 ( t ) y 2 ( t ) ,

where u 1 and u 2 are functions of t. Differentiating Eq. (6.14) gives us

Observe that it is pointless to search for solutions of the form y p = c 1 y 1 + c 2 y 2 where c 1 and c 2 are constants because for every choice of c 1 and c 2 , c 1 y 1 + c 2 y 2 is a solution to the corresponding homogeneous equation.

y p = u 1 y 1 + u 1 y 1 + u 2 y 2 + u 2 y 2 .

Assuming that

(6.15) y 1 u 1 + y 2 u 2 = 0

results in y p = u 1 y 1 + u 2 y 2 . Computing the second derivative then yields

y p = u 1 y 1 + u 1 y 1 + u 2 y 2 + u 2 y 2 .

Substituting y p , y p , and y p into Eq. (6.8) and using the facts that

u 1 ( y 1 + p y 1 + q y 1 ) = 0 and u 2 ( y 2 + p y 2 + q y 2 ) = 0

(because y 1 and y 2 are solutions to the corresponding homogeneous equation) results in

(6.16) d 2 y p d t 2 + p ( t ) d y p d t + q ( t ) y p = u 1 y 1 + u 1 y 1 + u 2 y 2 + u 2 y 2 + p ( t ) ( u 1 y 1 + u 2 y 2 ) + q ( t ) ( u 1 y 1 + u 2 y 2 ) = y 1 u 1 + y 2 u 2 = f ( t ) .

Observe that Eq. (6.15) and Eq. (6.16) form a system of two linear equations in the unknowns u 1 and u 2 :

(6.17) y 1 u 1 + y 2 u 2 = 0 y 1 u 1 + y 2 u 2 = f ( t ) .

Applying Cramer's Rule gives us

(6.18) u 1 = | 0 y 2 f ( t ) y 2 | | y 1 y 2 y 1 y 2 | = y 2 ( t ) f ( t ) W ( S ) and u 2 = | y 1 0 y 1 f ( t ) | | y 1 y 2 y 1 y 2 | = y 1 ( t ) f ( t ) W ( S ) ,

where W ( S ) is the Wronskian, W ( S ) = | y 1 y 2 y 1 y 2 | . After integrating to obtain u 1 and u 2 , we form y p , and then a general solution, y = y h + y p .

Example 6.15

Solve y + 9 y = sec 3 t , y ( 0 ) = 0 , y ( 0 ) = 0 , 0 t < π / 6 .

Solution

The corresponding homogeneous equation is y + 9 y = 0 with general solution y h = c 1 cos 3 t + c 2 sin 3 t . Then, a fundamental set of solutions is S = { cos 3 t , sin 3 t } and W ( S ) = 3 , as we see using Det, and Simplify.

fs = { Cos [ 3 t ] , Sin [ 3 t ] } ;

wm = { fs , D [ fs , t ] } ;

wm // MatrixForm

wd = Det [ wm ] // Simplify

( Cos [ 3 t ] Sin [ 3 t ] 3 Sin [ 3 t ] 3 Cos [ 3 t ] )

3

We use Eq. (6.18) to find u 1 = 1 9 ln cos 3 t and u 2 = 1 3 t .

u1 = Integrate [ Sin [ 3 t ] Sec [ 3 t ] / 3 , t ]

u2 = Integrate [ Cos [ 3 t ] Sec [ 3 t ] / 3 , t ]

1 9 Log [ Cos [ 3 t ] ]

t 3

It follows that a particular solution of the nonhomogeneous equation is y p = 1 9 cos 3 t ln cos 3 t + 1 3 t sin 3 t , and a general solution is y = y h + y p = c 1 cos 3 t + c 2 sin 3 t + 1 9 cos 3 t ln cos 3 t + 1 3 t sin 3 t .

Absolute value is not needed in the antiderivatives, because we are restricting the domain to 0 t &lt; π / 6 and cos t &gt; 0 on this interval.

yp = u1 Cos [ 3 t ] + u2 Sin [ 3 t ]

1 9 Cos [ 3 t ] Log [ Cos [ 3 t ] ] + 1 3 t Sin [ 3 t ]

Identical results are obtained using DSolve.

The negative sign in the output does not affect the result, because C[1] is arbitrary.

DSolve [ y " [ t ] + 9 y [ t ] = = Sec [ 3 t ] , y [ t ] , t ]

{ { y [ t ] C [ 1 ] Cos [ 3 t ] + C [ 2 ] Sin [ 3 t ] + 1 9 ( Cos [ 3 t ] Log [ Cos [ 3 t ] ] + 3 t Sin [ 3 t ] ) } }

Applying the initial conditions gives us c 1 = c 2 = 0 , so we conclude that the solution to the initial-value problem is y = 1 9 cos 3 t ln cos 3 t + 1 3 t sin 3 t .

sol = DSolve [ { y " [ t ] + 9 y [ t ] = = Sec [ 3 t ] , y [ 0 ] = = 0 , y [ 0 ] = = 0 } , y [ t ] , t ]

{ { y [ t ] 1 9 ( Cos [ 3 t ] Log [ Cos [ 3 t ] ] + 3 t Sin [ 3 t ] ) } }

We graph the solution with Plot in Fig. 6.18.

Figure 6.18

Figure 6.18. The domain of the solution is 0 ⩽t &lt;π/6 (University of Kentucky colors).

Plot [ y [ t ] /. sol , { t , 0 , Pi / 6 } ,

PlotStyle { { Thickness [ . 01 ] , CMYKColor [ 1 , . 75 , 0 , 0 ] } } ]  

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128241639000112

First Order Ordinary Differential Equations

Martha L. Abell , James P. Braselton , in Introductory Differential Equations (Fifth Edition), 2018

2.3.2.1 Outline of the Method of Undetermined Coefficients to Solve y + k y = q ( t ) Where k Is Constant and q ( t ) Is a Linear Combination of the Functions in (2.5)

1.

Solve the corresponding homogeneous equation, y + k y = 0 : y h ( t ) = C e k t .

2.

Determine the form of a particular solution y p ( t ) . (See Determining the Form of y p ( t ) next.)

3.

Determine the unknown coefficients in y p ( t ) by substituting y p ( t ) into the nonhomogeneous equation, y + k y = q ( t ) , and equating the coefficients of like terms.

4.

Form a general solution with y ( t ) = y h ( t ) + y p ( t ) = C e k t + y p .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128149485000021

Higher-Order Differential Equations

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

4.5.2 Higher-Order Nonhomogeneous Equations

In the same way as with second-order equations, we assume that a particular solution of the nth-order linear nonhomogeneous equation

y ( n ) + a n 1 ( t ) y ( n 1 ) + + a 2 ( t ) y + a 1 ( t ) y + a 0 ( t ) y = f ( t )

has the form y p = u 1 ( t ) y 1 + u 2 ( t ) y 2 + + u n ( t ) y n , where S = { y 1 , y 2 , , y n } is a fundamental set of solutions to the corresponding homogeneous equation

y ( n ) + a n 1 ( t ) y ( n 1 ) + + a 2 ( t ) y + a 1 ( t ) y + a 0 ( t ) y = 0 .

With the assumptions

(4.18) y p = y 1 u 1 + y 2 u 2 + + y n u n = 0 y p = y 1 u 1 + y 2 u 2 + + y n u n = 0 y p ( n 1 ) = y 1 ( n 2 ) u 1 + y 2 ( n 2 ) u 2 + + y n ( n 2 ) u n = 0

we obtain the equation

(4.19) y 1 ( n 1 ) u 1 + y 2 ( n 1 ) u 2 + + y n ( n 1 ) u n = f ( t ) .

Equations (4.18) and (4.19) form a system of n linear equations in the unknowns u 1′, u 2′, , u n ′. Applying Cramer's Rule,

(4.20) u i = W i ( S ) W ( S ) ,

where W(S) is given by equation (4.7),

W ( S ) = y 1 y 2 y n y 1 y 2 y n y 1 ( n 1 ) y 2 ( n 1 ) y n ( n 1 ) ,

and W i (S) is the determinant of the matrix obtained by replacing the ith column of

y 1 y 2 y n y 1 y 2 y n y 1 ( n 1 ) y 2 ( n 1 ) y n ( n 1 ) by 0 0 f ( t ) .

Example 4.5.2

Solve y ( 3 ) + 4 y = sec 2 t .

Solution

A general solution of the corresponding homogeneous equation is y h = c 1 + c 2 cos 2 t + c 3 sin 2 t ; a fundamental set is S = 1 , cos 2 t , sin 2 t with Wronskian W(S) = 8.

yh=DSolve[y"'[t]+4y′[t]==0, y[t], t]

{ { y [ t ] C [ 3 ] 1 2 C [ 2 ] Cos [ 2 t ] + 1 2 C [ 1 ] Sin [ 2 t ] } }

s={1, Cos[2t], Sin[2t]};Wronskian[s, t]

8

Using variation of parameters to find a particular solution of the nonhomogeneous equation, we let y 1 = 1, y 2 = cos 2 t , and y 3 = sin 2 t and assume that a particular solution has the form y p = u 1 y 1 + u 2 y 2 + u 3 y 3. Using the variation of parameters formula, we obtain

u 1 = 1 8 0 cos 2 t sin 2 t 0 2 sin 2 t 2 cos 2 t sec 2 t 4 cos 2 t 4 sin 2 t = 1 4 sec 2 t so u 1 = 1 8 ln | sec 2 t + tan 2 t | , u 2 = 1 8 1 0 sin 2 t 0 0 2 cos 2 t 0 sec 2 t 4 sin 2 t = 1 4 so u 2 = 1 4 t and u 3 = 1 8 1 cos 2 t 0 0 2 sin 2 t 0 0 4 cos 2 t sec 2 t = 1 4 tan 2 t so u 3 = 1 8 ln | cos 2 t | ,

where we use Det and Integrate to evaluate the determinants and integrals. In the case of u 1, the output given by Mathematica looks different than the result we obtained by hand but using properties of logarithms ( ln a / b = ln a ln b ) and trigonometric identities ( cos 2 x + sin 2 x = 1 , sin 2 x = 2 sin x cos x , cos 2 x sin 2 x = cos 2 x , and the reciprocal identities) shows us that

1 8 ln | cos t + sin t | ln | cos t + sin t | = 1 8 ln cos t + sin t cos t sin t = 1 8 ln cos t + sin t cos t sin t cos t + sin t cos t + sin t = 1 8 ln cos 2 t + 2 cos t sin t + sin 2 t cos 2 t sin 2 t = 1 8 ln 1 + sin 2 t cos 2 t = 1 8 ln 1 cos 2 t + sin 2 t cos 2 t = 1 8 ln sec 2 t + tan 2 t

so the results obtained by hand and with Mathematica are the same.

u1p=1/8Det[{{0, Cos[2t], Sin[2t]}, {0, −2Sin[2t], 2Cos[2t]}, {Sec[2t], −4Cos[2t], −4Sin[2t]}}]//Simplify

1 4 Sec [ 2 t ]

u1=Integrate[u1p, t]

1 4 ( 1 2 Log [ Cos [ t ] Sin [ t ] ] + 1 2 Log [ Cos [ t ] + Sin [ t ] ] )

u2p=1/8Det[{{1, 0, Sin[2t]}, {0, 0, 2Cos[2t]}, {0, Sec[2t], −4Sin[2t]}}]//Simplify

1 4

u2=Integrate[u2p, t]

t 4

u3p=1/8Det[{{1, Cos[2t], 0}, {0, −2Sin[2t], 0}, {0, −4Cos[2t], Sec[2t]}}]//Simplify

1 4 Tan [ 2 t ]

u3=Integrate[u3p, t]

1 8 Log [ Cos [ 2 t ] ]

Thus, a particular solution of the nonhomogeneous equation is

y p = 1 8 ln | sec 2 t + tan 2 t | 1 4 t cos 2 t + 1 8 ln | cos 2 t | sin 2 t

and a general solution is y = y h + y p . We verify the calculations using DSolve, which returns an equivalent solution.

gensol=DSolve[y"'[t]+4y′[t]==Sec[2t], y[t], t]//Simplify

{ { y [ t ] 1 8 ( 8 C [ 3 ] 8 C [ 2 ] Cos [ t ] 2 2 t Cos [ 2 t ] Log [ Cos [ t ] Sin [ t ] ] + Log [ Cos [ t ] + Sin [ t ] ] + 4 C [ 1 ] Sin [ 2 t ] + Log [ Cos [ 2 t ] ] Sin [ 2 t ] ) } }

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000048

Higher Order Equations

Martha L. Abell , James P. Braselton , in Introductory Differential Equations (Fourth Edition), 2014

Undetermined Coefficients

In the special case that the corresponding homogeneous equation has constant coefficients and the forcing function is a linear combination of functions of the form

(4.23) 1 , t , , t 2 , , e α t , t e α t , t 2 e α t , , cos β t , sin β t , t cos β t , t sin β t , t 2 cos β t , t 2 sin β t , or e α t cos β t , e α t sin β t , t e α t cos β t , t e α t sin β t , t 2 e α t cos β t , t 2 e α t sin β t , ,

then the method of undetermined coefficients can be used to determine the form of a particular solution of the nonhomogeneous equation in the same way as that discussed in Section 4.2.

Example 4.6.2

Solve y (5) + 4y‴ = 48t − 6 − 10et .

Solution

The corresponding homogeneous equation is y (5) + 4y‴ = 0, which has characteristic equation r 5 + 4r 3 = r 3(r 2 + 4) = 0 so r 1,2,3 = 0 has multiplicity three and r 4,5 = ±2i each have multiplicity one. Thus, a fundamental set for the corresponding homogeneous equation is S = 1 , t , t 2 , cos 2 t , sin 2 t and a general solution of the corresponding homogeneous equation is y h = c 1 + c 2 t + c 3 t 2 + c 4 cos 2 t + c 5 sin 2 t .

For the forcing function, f(t) = 48t − 6 − 10et we have two associated sets:

F 1 = e t and F 2 = t , 1 ,

corresponding to the terms − 10et and 48t − 6, respectively, in the forcing function. Note that no element of F 1 is a solution of the corresponding homogeneous equation. On the other hand, functions in F 2 are solutions to the corresponding homogeneous equation. So, following the outline in Section 4.2, we multiply F 2 by t n where n is the smallest positive integer so that no function in t n F 2 is a solution of the corresponding homogeneous equation. In this case, we multiply F 2 by t 3 to obtain t 3 F 2 = t 4 , t 3 .

Thus, we assume that a particular solution to the nonhomogeneous equation has the form y p = At 4 + Bt 3 + Cet . Differentiating we have,

y p = 4 A t 3 + 3 B t 2 C e t , y p = 1 2 A t 2 + 6 B t + C e t , y p = 2 4 A + C e t and y p ( 4 ) = C e t .

Substituting into the nonhomogeneous equation, simplifying the result, and equating coefficients gives us

2 4 B + 9 6 A t 5 C e t = 4 8 t 6 1 0 e t ,

so 24B = −6, 96A = 48, and − 5C = −10. Thus, A = 1/2, B = −1/4, and C = 2 so a particular solution of the nonhomogeneous equation is y p = 1 2 t 4 1 4 t 3 + 2 e t .

A general solution of the nonhomogeneous equation is then given by

y = y h + y p = c 1 + c 2 t + c 3 t 2 + c 4 cos 2 t + c 5 sin 2 t + 1 2 t 4 1 4 t 3 + 2 e t .

What is the form of a particular solution of y ( 5 ) + 4 y = t cos 2 t ?

To solve an IVP, first determine a general solution and then use the initial conditions to solve for the unknown constants in the general solution.

Example 4.6.3

Solve the IVP

4 y + 4 y + 6 5 y = e t 3 2 8 6 3 cos 2 t + 5 7 7 2 7 sin 2 t ,

y(0) = 4/3, y′(0) = −5/2, y″(0) = −173/12.

Solution

The corresponding homogeneous equation is 4 y‴ + 4y″ + 65y′ = 0 with characteristic equation 4r 3 + 4r 2 + 65r = r(4r 2 + 4r + 65) = 0 so r 1 = 0 and using the quadratic formula to solve 4r 2 + 4r + 65 = 0 gives us r 2 , 3 = 1 2 ± 4 i . Thus, a fundamental set of solutions for the corresponding homogeneous equation is S = { 1 , e t 2 cos 4 t , e t 2 sin 4 t } and a general solution of the corresponding homogeneous equation is y h = c 1 + e t 2 ( c 2 cos 4 t + c 3 sin 4 t ) .

The associated set of functions for the forcing function

g ( t ) = e t 3 2 8 6 3 cos 2 t + 5 7 7 2 7 sin 2 t

is F = e t 3 cos 2 t , e t 3 sin 2 t . No function in F is a solution of the corresponding homogeneous equation so we assume that a particular solution of the nonhomogeneous equation has the form y p = A e t 3 cos 2 t + B e t 3 sin 2 t , where a and b are constants to be determined. Differentiating y p three times gives us

y p = 1 3 A + 2 B e t 3 cos 2 t + 2 A 1 3 B e t 3 sin 2 t , y p = 3 5 9 A 4 3 B e t 3 cos 2 t + 4 3 A 3 5 9 e t 3 sin 2 t , and y p = 1 0 7 2 7 A 2 2 3 B e t 3 cos 2 t + 2 2 3 A + 1 0 7 2 7 B e t 3 sin 2 t .

Substituting y p and its derivatives into the nonhomogeneous equation and simplifying the result gives us

When algebra becomes unusually cumbersome, we usually use a computer algebra system to assist in checking our calculations.

4 y p + 4 y p + 6 5 y p = 5 7 7 2 7 A + 2 8 6 3 B e t 3 cos 2 t + 2 8 6 3 A 5 7 7 2 7 B e t 3 sin 2 t = e t 3 2 8 6 3 cos 2 t + 5 7 7 2 7 sin 2 t .

Equating coefficients gives us the system

5 7 7 2 7 A + 2 8 6 3 B = 2 8 6 3 2 8 6 3 A 5 7 7 2 7 B = 5 7 7 2 7 ,

which has solution A = 0 and B = −1. Thus, y p = e t 3 sin 2 t and a general solution of the nonhomogeneous equation is y = y h + y p = c 1 + e t 2 ( c 2 cos 4 t + c 3 sin 4 t e t 3 sin 2 t (see Figure 4.11(a)).

Figure 4.11. (a) Various solutions of the nonhomogeneous equation. (b) The solution of the nonhomogeneous equation that satisfies the initial conditions y(0) = 4/3, y (0) = −5/2, and y (0) = −173/12.

To solve the IVP, we first differentiate Y twice resulting in

y = 1 2 c 2 + 4 c 3 e t 2 cos 4 t + 4 c 2 1 2 c 3 e t 2 sin 4 t 2 e t 3 cos 2 t + 1 3 e t 3 sin 2 t and y = 6 3 4 c 2 + 4 c 3 e t 2 cos 4 t + 4 c 2 6 3 4 c 3 e t 2 sin 4 t + 4 3 e t 3 cos 2 t + 3 5 9 e t 3 sin 2 t .

Evaluating at t = 0 gives us the system of equations

y ( 0 ) = c 1 + c 2 = 4 3 y ( 0 ) = 2 1 2 c 2 + 4 c 3 = 5 2 y ( 0 ) = 2 2 3 + 1 9 1 8 c 2 6 1 c 3 = 1 7 3 1 2 ,

which has solution c 1 = 1/3, c 2 = 1, and c 3 = 0 so the solution to the IVP is y = 1 3 + e t 2 cos 4 t e t 3 sin 2 t (see Figure 4.11(b)).
What is the form of a particular solution of 4y‴ + 4y″ + 65y′ = t2 + te-t/3 cos2t?

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780124172197000041

Laplace Transform Methods

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

Abstract

In previous chapters we have investigated solving the nth-order linear equation

a n ( t ) y ( n ) + a n 1 ( t ) y ( n 1 ) + + a 2 ( t ) y + a 1 ( t ) y + a 0 ( t ) y = f ( t )

for y. We have seen that if the coefficients a i (t ) are numbers, we can find a general solution of the equation by first solving the characteristic equation of the corresponding homogeneous equation, forming a general solution of the corresponding homogeneous equation, and then finding a particular solution to the nonhomogeneous equation. If the coefficients a i (t) are not constants, we have learned that solving above equation may be substantially more difficult that in other cases, such as when the functions a i (t) are constants. For example, when above equation is a Cauchy-Euler equation, techniques used to solve the case when above equation has constant coefficients can be used to solve the equation. In other situations, we might be able to use a series to find a solution of the equation. Regardless, in all these situations the function f(t) has (typically) been a smooth function. If f(t) is not a smooth function, such as when f(t) is a piecewise-defined or periodic function, solving above equation can be substantially more difficult. In this chapter, we discuss a technique that transforms above equation into an algebraic equation that can sometimes be solved so that a solution to the differential equation can be obtained.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000085

Green's Functions

George B. Arfken , ... Frank E. Harris , in Mathematical Methods for Physicists (Seventh Edition), 2013

Form of Green's Function

The properties we have identified for G are sufficient to enable its more complete identification, given a Hermitian operator and its boundary conditions. We continue with the study of problems on an interval (a, b) with one homogeneous boundary condition at each endpoint of the interval.

Given a value of t, it is necessary for x in the range ax < t that G (x, t) have an x dependence y 1(x) that is a solution to the homogeneous equation = 0 and that also satisfies the boundary condition at x = a. The most general G (x, t) satisfying these conditions must have the form

(10.16) G ( x , t ) = y 1 ( x ) h 1 ( t ) , ( x < t ) ,

where h 1(t) is presently unknown. Conversely, in the range t < xb, it is necessary that G (x, t) have the form

(10.17) G ( x , t ) = y 2 ( x ) h 2 ( t ) , ( x > t ) ,

where y 2 is a solution of = 0 that satisfies the boundary condition at x = b. The symmetry condition, Eq. (10.15), permits Eqs. (10.16) and (10.17) to be consistent only if h 2 = A y 1 and h 1 = A y 2 , with A a constant that is still to be determined. Assuming that y 1 and y 2 can be chosen to be real, we are led to the conclusion that

(10.18) G ( x , t ) = A y 1 ( x ) y 2 ( t ) , x < t , A y 2 ( x ) y 1 ( t ) , x > t ,

where y i = 0 , with y 1 satisfying the boundary condition at x = a and y 2 satisfying that at x = b. The value of A in Eq. (10.18) depends, of course, on the scale at which the yi have been specified, and must be set to a value that is consistent with Eq. (10.9). As applied here, that condition reduces to

A y 2 ( t ) y 1 ( t ) y 1 ( t ) y 2 ( t ) = 1 p ( t ) ,

equivalent to

(10.19) A = p ( t ) [ y 2 ( t ) y 1 ( t ) y 1 ( t ) y 2 ( t ) ] 1 .

Despite its appearance, A does not depend on t. The expression involving the yi is their Wronskian, and it has a value proportional to 1/p(t). See Exercise 7.6.11.

It is instructive to verify that the form for G(x, t) given by Eq. (10.18) causes Eq. (10.7) to generate the desired solution to the ODE y = f . To this end, we obtain an explicit form for y(x):

(10.20) y ( x ) = A y 2 ( x ) a x y 1 ( t ) f ( t ) d t + A y 1 ( x ) x b y 2 ( t ) f ( t ) d t .

From Eq. (10.20) it is easy to verify that the boundary conditions on y(x) are satisfied; if x = a the first of the two integrals vanishes, and the second is proportional to y 1; corresponding remarks apply at x = b.

It remains to show that Eq. (10.20) yields y = f . Differentiating with respect to x, we first have

(10.21) y ( x ) = A y 2 ( x ) a x y 1 ( t ) f ( t ) d t + A y 2 ( x ) y 1 ( x ) f ( x ) + A y 1 ( x ) x b y 2 ( t ) f ( t ) d t A y 1 ( x ) y 2 ( x ) f ( x ) = A y 2 ( x ) a x y 1 ( t ) f ( t ) d t + A y 1 ( x ) x b y 2 ( t ) f ( t ) d t .

Proceeding to (py′)′:

(10.22) [ p ( x ) y ( x ) ] = A [ p ( x ) y 2 ( x ) ] a x y 1 ( t ) f ( t ) d t + A [ p ( x ) y 2 ( x ) ] y 1 ( x ) f ( x ) + A [ p ( x ) y 1 ( x ) ] x b y 2 ( t ) f ( t ) d t A [ p ( x ) y 1 ( x ) ] y 2 ( x ) f ( x ) .

Combining Eq. (10.22) and q(x) times Eq. (10.20), many terms drop because y 1 = y 2 = 0 , leaving

(10.23) y ( x ) = A p ( x ) y 2 ( x ) y 1 ( x ) y 1 ( x ) y 2 ( x ) f ( x ) = f ( x ) ,

where the final simplification took place using Eq. (10.19).

Example 10.1.1

Simple Second-Order ODE

Consider the ODE

y = f ( x ) ,

with boundary conditions y(0) = y (1) = 0. The corresponding homogeneous equation y″ = 0 has general solution y 0 = c 0 + c 1 x; from these we construct the solution y 1 = x that satisfies y 1(0) = 0 and the solution y 2 = 1 − x, satisfying y 2(1) = 0. For this ODE, the coefficient p(x) = −1, y 1 ( x ) = 1 , y 2 ( x ) = 1 , and the constant A in the Green's function is

A = [ ( 1 ) , [ ( 1 ) ( x ) ( 1 ) ( 1 x ) ] ] 1 = 1.

Our Green's function is therefore

G ( x , t ) = x ( 1 t ) , 0 x < t , t ( 1 x ) , t < x 1 .

Assuming we can perform the integral, we can now solve this ODE with boundary conditions for any function f (x). For example, if f (x) = sin πx, our solution would be

y ( x ) = 0 1 G ( x , t ) sin π t d t = ( 1 x ) 0 x t sin π t d t + x x 1 ( 1 t ) sin π t d t = 1 π 2 sin π x .

The correctness of this result is easily checked.

One advantage of the Green's function formalism is that we do not need to repeat most of our work if we change the function f (x). If we now take f (x) = cos πx, we get

y ( x ) = 1 π 2 ( 2 x 1 + cos π x ) .

Note that our solution takes full account of the boundary conditions.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123846549000104

First-Order Ordinary Differential Equations

Martha L. Abell , James P. Braselton , in Differential Equations with Mathematica (Fourth Edition), 2016

2.5.2 Variation of Parameters and the Method of Undetermined Coefficients

Observe that equation (2.4) is separable:

d y d x + p ( x ) y = 0 1 y d y = p ( x ) d x ln y = p ( x ) d x + C y = C e p ( x ) d x .

Notice that any constant multiple of a solution to a linear homogeneous equation is also a solution. Now suppose that y is any solution of equation (2.3) and y p is a particular solution of equation (2.3). Then,

A particular solution is a specific function that is a solution to the equation that does not contain any arbitrary constants.

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

Thus, yy p is a solution to the corresponding homogeneous equation of equation ( 2.3). Hence,

y y p = C e p ( x ) d x y = C e p ( x ) d x + y p y = y h + y p ,

where y h = C e p ( x ) d x . That is, a general solution of equation (2.3) is

y = y h + y p ,

where y p is a particular solution to the nonhomogeneous equation and y h is a general solution to the corresponding homogeneous equation. Thus, to solve equation (2.3), we need to first find a general solution to the corresponding homogeneous equation, y h , which we can accomplish through separation of variables, and then find a particular solution, y p , to the nonhomogeneous equation.

If y h is a solution to the corresponding homogeneous equation of equation (2.3) then for any constant C, Cy h is also a solution to the corresponding homogeneous equation. Hence, it is impossible to find a particular solution to equation (2.3) of this form. Instead, we search for a particular solution of the form y p = u(x)y h , where u(x) is not a constant function. Assuming that a particular solution, y p , to equation (2.3) has the form y p = u(x)y h , differentiating gives us

y p = u y h + u y h

and substituting into equation (2.3) results in

y p + p ( x ) y p = u y h + u y h + p ( x ) u y h = q ( x ) .

Because u y h + p ( x ) u y h = u y h + p ( x ) y h = u 0 = 0 , we obtain

y h is a solution to the corresponding homogeneous equation so y h ′ + p(x)y h = 0.

u y h = q ( x ) u = 1 y h q ( x ) u = e p ( x ) d x q ( x ) u = e p ( x ) d x q ( x ) d t

so

y p = u ( x ) y h = C e p ( x ) d x e p ( x ) d x q ( x ) d x .

Because we include an arbitrary constant of integration when evaluating e p ( x ) d x q ( x ) d x , it follows that we can write a general solution of equation (2.3) as

(2.6) y = e p ( x ) d x e p ( x ) d x q ( t ) d t .

Example 2.5.5 (Exponential Growth)

Let y = y(t) denote the size of a population at time t. If y grows at a rate proportional to the amount present, y satisfies

Exponential growth is discussed in more detail in Section 3.2.1.

(2.7) d y d t = α y ,

where α is the growth constant. If y(0) = y 0, using equation (2.6) results in y = y 0 e αt . We use DSolve to confirm this result.

Clear[t, y]DSolve[{y′[t]==αy[t], y[0]==y0}, y[t], t]

y [ t ] e t α y0

d y / d t = k y y s models Newton's Law of Cooling: the rate at which the temperature, y(t), changes in a heating/cooling body is proportional to the difference between the temperature of the body and the constant temperature, y s , of the surroundings. Newton's Law of Cooling is discussed in more detail in Section 3.3.

Example 2.5.6

Solve each of the following equations: (a) d y / d t = k y y s , y(0) = y 0, k and y s constant and (b) y′− 2ty = t.

Solution

(a) By hand, we rewrite the equation and obtain

d y d t k y = k y s .

A general solution of the corresponding homogeneous equation

d y d t k y = 0

is y h = e kt . Because k and − ky s are constants, we suppose that a particular solution of the nonhomogeneous equation, y p , has the form y p = A, where A is a constant.

Assuming that y p = A, we have y p ′ = 0 and substitution into the nonhomogeneous equation gives us

This will turn out to be a lucky guess. If there is not a solution of this form, we would not find one of this form.

d y p d t k y p = K A = k y s so A = y s .

Thus, a general solution is y = y h + y p = Ce kt + y s . Applying the initial condition y(0) = y 0 results in y = y s + (y 0y s )e kt .

We obtain the same result with DSolve. We graph the solution satisfying y(0) = 75 assuming that k = −1/2 and y s = 300 in Figure 2-23. Notice that y(t) → y s as t .

Figure 2-23. The temperature of the body approaches the temperature of its surroundings

sola=DSolve[{y′[t]==k(y[t]−ys), y[0]==y0}, y[t], t] y [ t ] e k t y0 + ys e k t ys

tp=sola[[1, 1, 2]]/.{k→−1/2, ys→300, y0→75};p1=Plot[tp, {t, 0, 10}, AxesLabel→{t, y}, AxesStyle→Black]

(b) The equation is in standard form and we identify p(t) = −2t. Then, the integrating factor is μ ( t ) = e p ( t ) d t = e t 2 . Multiplying the equation by the integrating factor, μ(t), results in

e t 2 ( y 2 t y ) = t e t 2 or d d t y e t 2 = t e t 2 .

Integrating gives us

y e t 2 = 1 2 e t 2 + C or y = 1 2 + C e t 2 .

We confirm the result with DSolve.

solb=DSolve[y′[t]−2ty[t]==t, y[t], t]

y [ t ] 1 2 + e t 2 C [ 1 ]

Application: Antibiotic Production

When you are injured or sick, your doctor may prescribe antibiotics to prevent or cure infections. In the journal article "Changes in the Protein Profile of Streptomyces Griseus during a Cycloheximide Fermentation" we see that production of the antibiotic cycloheximide by Streptomyces is typical of antibiotic production. During the production of cycloheximide, the mass of Streptomyces grows relatively quickly and produces little cycloheximide. After approximately 24 hours, the mass of Streptomyces remains relatively constant and cycloheximide accumulates. However, once the level of cycloheximide reaches a certain level, extracellular cycloheximide is degraded (feedback inhibited). One approach to alleviating this problem to maximize cycloheximide production is to continuously remove extracellular cycloheximide. The rate of growth of Streptomyces can be described by the separable equation

Source: Kevin H. Dykstra and Henry Y. Wang, "Changes in the Protein Profile of Streptomyces Griseus during a Cycloheximide Fermentation," Biochemical Engineering V, Annals of the New York Academy of Sciences, Volume 56, New York Academy of Sciences (1987), pp. 511–522.

d X d t = μ max 1 1 X max X X ,

where X represents the mass concentration in g/L, μ max is the maximum specific growth rate, and X max represents the maximum mass concentration. We now solve the initial-value problem d X / d t = μ max 1 1 X max X X X ( 0 ) = 1 with DSolve, naming the result sol1.

Note that this equation can be converted to a linear equation with the substitution y = X −1.

Clear[x]sol1=DSolve[{x′[t]==μ(1−x[t]/xmax)x[t], x[0]==1}, x[t], t]

x [ t ] e t μ xmax 1 + e t μ + xmax

Experimental results have shown that μ max = 0.3 h r 1 and X max = 10 g / L . For these values, we use Plot to graph X(t) on the interval [0, 24] in Figure 2-24 (a). Then, we use Table and TableForm to determine the mass concentration at the end of 4, 8, 12, 16, 20, and 24 hours.

Figure 2-24. (a) Plot of the mass concentration, x(t). (b) Accumulation of the antibiotic

μ=0.3;xmax=10.;

p1=Plot[x[t]/.sol1, {t, 0, 24}, PlotStyle→Black, AxesStyle→Black, AxesLabel→{t, x}, PlotLabel→"(a)"];

TableForm[Table[{t, sol1[[1, 1, 2]]}, {t, 4, 24, 4.}]]

4 . 2.69487 8 . 5.50521 12 . 8.02624 16 . 9.3104 20 . 9.78178 24 . 9.93326

The rate of accumulation of cycloheximide is the difference between the rate of synthesis and the rate of degradation:

d P d t = R s R d .

It is known that R d = K d P, where K d = 5 × 10−3 h−1, so dP/dt = R s R d is equivalent to dP/dt = R s K d P. Furthermore,

R s = Q p o E X 1 + P K l 1 ,

where Q po represents the specific enzyme activity with value Q po ≈ 0.6   g CH/g,protein ⋅h and K l represents the inhibition constant. E represents the intracellular concentration of an enzyme which we will assume is constant. For large values of K l and t, X(t) ≈ 10 and 1 + P / K l 1 1 . Thus, R s ≈ 10Q po E so

d P d t = 10 Q p o E K d P .

After defining K d and Q po , we solve the initial-value problem d P / d t = 10 Q p o E K d P p ( 24 ) = 0 and then graph 1 E P ( t ) on the interval [0, 24] in Figure 2-24 (b).

Clear[p]kd=5/1000;qp0=0.6;sol2=DSolve[{p′[t]==10qp0−kdp[t], p[24]==0}, p[t], t]

p [ t ] 1200 . e 0.005 t 1.1275 + 1 . e 0.005 t

p2=Plot[p[t]/.sol2, {t, 24, 1000}, PlotStyle→Black, AxesStyle→Black, AxesLabel→{t, p}, AxesOrigin→{0, 0}, PlotLabel→"(b)"];

Show[GraphicsRow[{p1, p2}]]

From the graph, we see that the total accumulation of the antibiotic approaches a limiting value, which in this case is 1200.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128047767000024

Ordinary differential equations

Brent J. Lewis , ... Andrew A. Prudil , in Advanced Mathematics for Engineering Students, 2022

Homogeneous linear equations

An ordinary differential equation of nth order is of the form

(2.67) F ( x , y , y , . . . y ( n ) ) = 0 ,

where the nth derivative y ( n ) = d n y d x n . The equation is linear if it can be written as

(2.68)

The corresponding homogeneous equation written in "standard" form is

(2.69)

Eq. (2.69) has the important property that a linear combination of solutions is again a solution (superposition principle or linearity principle). A basis of solutions of Eq. (2.69) consists of n linearly independent solutions of y 1 , . . . , y n . The corresponding general solution is a linear combination:

(2.70)

From Eq. (2.70) one obtains a particular solution by choosing n numbers for c 1 , . . . c n by imposing n initial conditions:

(2.71)

Together Eq. (2.69) and Eq. (2.71) constitute an initial value problem for Eq. (2.69). If p 0 , . . . , p n 1 are continuous on some open interval I and x 0 is in I, then Eq. (2.69) has a general solution on I, and Eq. (2.69) and Eq. (2.71) have a unique solution on I (which is a particular solution).

Similarly, as in Section 2.2.4, the basis functions y 1 , . . . , y n are linearly independent on the interval I if and only if the Wronskian (extended to nth order), as defined by

(2.72)

is different from zero on the interval.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128236819000101

Applications of Higher Order Differential Equations

Martha L. Abell , James P. Braselton , in Introductory Differential Equations (Fourth Edition), 2014

L-R-C Circuits

Second-order nonhomogeneous linear ODEs arise in the study of electrical circuits after the application of Kirchhoff's law. Suppose that I(t) is the current in the L-R-C series electrical circuit (shown in Figure 5.14) where L, r, and c represent the inductance, resistance, and capacitance of the circuit, respectively.

Figure 5.14. An L-R-C circuit.

The voltage drops across the circuit elements in Table 5.2 have been obtained from experimental data, where Q is the charge of the capacitor and dQ/dt = I.

Table 5.2. Voltage Drops Across an L-R-C Circuit

Circuit Element Voltage Drop
Inductor

L d I d t

Resistor RI
Capacitor

1 C Q

Our goal is to model this physical situation with an IVP so that we can determine the current and charge in the circuit. For convenience, the terminology used in this section is summarized in Table 5.3.

Table 5.3. Terminology Used in Section 5.4

Electrical Quantities Units
Inductance (L) Henrys (H)
Resistance (r) Ohms (ω)
Capacitance (c) Farads (F)
Charge (Q) Coulombs (C)
Current (I) Amperes (A)

The physical principle needed to derive the differential equation that models the L-R-C series circuit is Kirchhoff's law.

The German physicist Gustav Robert Kirchhoff (1824-1887) worked in spectrum analysis, optics, and electricity.

Kirchhoff's law The sum of the voltage drops across the circuit elements is equivalent to the voltage E(t) impressed on the circuit.

Applying Kirchhoff's law with the voltage drops in Table 5.2 yields the differential equation L d I d t + R I + C 1 Q = E ( t ) . Using the fact that dQ/dt = I, we also have d2 Q/dt 2 = dI/dt. Therefore, the equation becomes L d 2 Q d t 2 + R d Q d t + C 1 Q = E ( t ) , which can be solved by the method of undetermined coefficients or the method of variation of parameters. If the initial charge and current are Q(0) = Q 0 and I(0) = Q'(0) = I 0, we solve the IVP

(5.8) L d 2 Q d t 2 + R d Q d t + 1 C Q = E ( t ) , Q ( 0 ) = Q 0 , I ( 0 ) = d Q d t ( 0 ) = I 0

for the charge Q(t). The solution is differentiated to find the current I(t).

Example 5.4.1

Consider the L-R-C circuit with L = 1 H, R = 40ω, C = 1/4000 F, and E(t) = 24 V. Determine the current in this circuit if there is zero initial current and zero initial charge.

Solution

Using the indicated values, the IVP that we must solve is

d 2 Q d t 2 + 4 0 d Q d t + 4 0 0 0 Q = 2 4 , Q ( 0 ) = Q 0 , I ( 0 ) = d Q d t ( 0 ) = I 0

The characteristic equation of the corresponding homogeneous equation is r 2 + 40r + 4000 = 0 with roots r 1,2 = −20 ± 60i, so a general solution of the corresponding homogeneous equation is Q h ( t ) = e 2 0 t ( c 1 cos 6 0 t + c 2 sin 6 0 t ) . Because the voltage is the constant function E(t) = 24, we assume that the particular solution has the form Q p(t) = A. Substitution into d2 Q/dt 2 + 40 dQ/dt + 4000Q = 24 yields 4000A = 24 or A = 3/500. Therefore, a general solution of the nonhomogeneous equation is

Q ( t ) = Q h ( t ) + Q p ( t ) = e 2 0 t ( c 1 cos 6 0 t + c 2 sin 6 0 t ) + 3 5 0 0

with derivative

d Q d t ( t ) = 2 0 e 2 0 t ( c 1 cos 6 0 t + c 2 sin 6 0 t ) + 6 0 e 2 0 t ( c 1 sin 6 0 t + c 2 cos 6 0 t ) .

After application of the initial conditions, we have Q(0) = c 1 + 3/500 = 0 and dQ/dt(0) = −20c 1 + 60c 2 = 0. Therefore, c 1 = −3/500 and c 2 = −1/500, so the charge in the circuit is

Q ( t ) = e 2 0 t 3 5 0 0 cos 6 0 t 1 5 0 0 sin 6 0 t + 3 5 0 0 ,

and the current is given by

d Q d t ( t ) = 2 0 e 2 0 t 3 5 0 0 cos 6 0 t 1 5 0 0 sin 6 0 t + 6 0 e 2 0 t 3 5 0 0 sin 6 0 t 1 5 0 0 cos 6 0 t = 2 5 e 2 0 t sin 6 0 t .

These results indicate that in time the charge approaches the constant value of 3/500, which is known as the steady-state charge. Also, because of the exponential term, the current approaches zero as t increases. We show the graphs of Q(t) and I(t) in Figure 5.15(a) and (b) to verify these observations.

In Example 5.4.1 , how is the charge Q(t) affected if E(t) = 48 V? What happens to Q(t) if R = 4 0 1 0 ω?

Figure 5.15. Plots of the charge and current in an L-R-C circuit.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780124172197000053

Linear Operator Equations

William L. Dunn , J. Kenneth Shultis , in Exploring Monte Carlo Methods, 2012

Solution for a constant source term

To obtain the form of the solution to

(8.93) 2 u ( r , θ ) r 2 + 1 r u ( r , θ ) r + 1 r 2 2 u ( r , θ ) θ 2 β 2 u ( r , θ ) = q o ,

observe that up = q o2 is a particular solution. Addition of up to the most general physically realistic solution of the corresponding homogeneous equation gives the general solution of Eq. (8.93), namely,

(8.94) u ( r , θ ) R Θ = a 0 I 0 ( β r ) + = 1 I ( β r ) [ a cos θ + b sin θ ] + q o / β 2 .

To see how u(r) at the center of a circle is related to the value on the circumference, integrate Eq. (8.94) over θ to obtain

(8.95) 0 2 π u ( r , θ ) d θ = 2 π a 0 I 0 ( β r ) + 2 π q o / β 2 .

The value at the circle center u(0, θ) = u(r0 ) is independent of θ and, from Eq. (8.94), has the value

(8.96) u ( r 0 ) = a 0 + q o / β 2 .

Combining these two results to eliminate a0 gives the important result

(8.97) u ( r 0 ) = 1 I 0 ( β r ) 1 2 π 0 2 π u ( r , θ ) d θ + Q ( r ) ,

where

(8.98) Q ( r ) q o β 2 ( 1 1 I 0 ( β r ) ) .

From this result it is seen that u at the center of a circle of radius r equals Q(r) plus the average of u(r, θ)/I0r) on the circumference. Thus, the previous random walking-on-circles procedure is altered to score the particle's weight times Q(r) at each step, in addition to scoring the surface contribution at the end of the walk. After N random walks the value of u(r0) is then estimated as

(8.99) u ( r 0 ) 1 N i = 1 N { j = 1 n i W j ( i ) Q ( r j ) + W n i ( i ) ψ ( r n i * ) } .

In Example 8.4, this walking-on-circles procedure is used to solve a particular diffusion problem.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780444515759000087