Mathcad - lesson8.mcd

Unit III Lesson 8
Boundary Value Problems
The preceding lesson was about dynamic systems that were modeled by initial value problems.
Often, however, the mathematical description of a physical system demands that we solve a DE
that has ‘boundary conditions’. That means we know values of the function, one or more of its
derivatives, or even some linear combination of the function and its derivatives at two or more
points. Several such examples follow.
Deflection of a Beam
Many structures built by humans are constructed from girders or beams which deflect or distort
(Shall we say bend?) under their own weight or due to some external force. This deflection is
governed by a linear 4th order DE. The diagram here shows the deflection curve of a beam:
Suppose the x axis coincides with the axis of symmetry and the deflection is positive in a
downward direction. In the theory of elasticity, it is known that the bending moment M ( x) at a
point x along the axis of the beam is related to the load per unit length w( x) by the DE:
2
d M
2
= w( x)
dx
In addition, the bending moment M ( x) = E⋅ I⋅ k where E and I are constants and k is the curvature
of the elastic curve. E involves the elasticity of the beam's material (Steel is more rigid than lead,
for example) and I is the moment of inertia of a cross section of the beam about the ‘neutral
axis’. The product E⋅ I is called the flexural rigidity of the beam.
You may recall from calculus class that curvature k is given by the formula k =
y''
3
(1 + y'2)
2
When the deflection is small y’approaches 0 so the denominator approaches 1, so k = y'' is
used as the equation for curvature. Then M ( x) = E⋅ I⋅ k can be written M ( x) = E⋅ I⋅ y''.
2
Differentiating this twice gives the differential equation
d M
2
dx
4
= E⋅ I⋅
d y
4
dx
2
and since w( x) =
d M
2
dx
4
we have the 4th order differential equation w( x) = E⋅ I⋅
d y
4
dx
Boundary conditions associated with the equation depend on how the ends of the beam are
supported. A cantilever beam is fixed at one end and free at the other as in the diagram.
For an embedded cantilever beam, the deflection y( x) must satisfy two boundary conditions at
the embedded end and two at the free end:
y( 0 ) = 0 since there is no deflection at the embedded end.
y'( 0 ) = 0 since the deflection curve at the embedded end is tangent to the x axis
y''( L) = 0 since the bending moment is zero at the end
y'''( L) = 0 since the shear force at the end is zero.
If the beam is simply hinged on the fixed end, then y( 0 ) = 0 and y''( 0 ) = 0
If the beam is embedded at both ends, there is no vertical deflection and the line of deflection is
horizontal at the endpoints, which gives boundary conditions:
y( 0 ) = 0, y'( 0 ) = 0, y( L) = 0, y'( L) = 0. A picture of the situation is given here:
ex 1 A beam of length L is embedded at both ends, and the constant load w is uniform along its
0
length. (I.e. w( x) = w for all x in the interval ( 0 , L)). Find the function representing the deflection
0
of the beam.
4
Solving the differentiual equation w( x) = E⋅ I⋅
d y
just require 4 integrations after dividing both
4
dt
4
w
d y
sides by E⋅ I and replacing w( x) with w to get
0
=
4
dt
d
y=
dx
d
2
:
0
E⋅ I
⋅x + c
1
2
0 x
w
2
dx
E⋅ I
w
3
3
0
y=
⋅
+ c ⋅x + c
1
2
E⋅ I 2
3
w
0 x
2
d
y=
⋅
+ c ⋅x + c ⋅x + c
1
2
3
E
⋅
I
6
dx
4
w
0 x
3
2
y=
⋅
+ c x + c x + c x+ c
1
2
3
4
E⋅ I 24
(Noting that the constants c and c absorbed the fractions)
1
2
Boundary conditions y( 0 ) = 0 and y'( 0 ) = 0 force c = 0 and c = 0 so the solution has become
3
w
4
4
w
0 x
0
3
2
3
2
⋅
+ c ⋅ x + c ⋅ x with y' =
⋅ x + 3 ⋅ c ⋅ x + 2 ⋅ c ⋅ x noting that though c and c are
1
2
1
2
1
2
E⋅ I 24
6 ⋅ E⋅ I
indeed arbitrary but the relationship between the function and its derivative must be maintained.
Boundary conditions y( L) = 0 and y'( L) = 0 gives a system of equations
4
w
0 L
3
2
0=
⋅
+ c ⋅L + c ⋅L
1
2
E⋅ I 24
w
0
3
2
0=
⋅ L + 3⋅ c ⋅ L + 2 ⋅ c ⋅ L
1
2
6 ⋅ E⋅ I
or rewritten
4
w
0 L
3
2
−
⋅
= c ⋅L + c ⋅L
1
2
E⋅ I 24
w
0
3
2
−
⋅ L = 3⋅ c ⋅ L + 2 ⋅ c ⋅ L
1
2
6 ⋅ E⋅ I
y=
Which can be solved for c and c using Cramer's rule (unit III lesson 5):
1
2
 w0 L4 2 
− ⋅
L 
 E⋅ I 24

 w

5
w0⋅ L
 − 0 ⋅ L3 2L 
−w ⋅ L
 6 ⋅ E⋅ I
 = 12⋅ E⋅ I = 0
c =
1
4
12⋅ E⋅ I
 L3 L2 
−L


 2

 3L 2L 
4 
w
 3
0
 L − ⋅L 
E⋅ I 24 



6
w
− w0⋅ L
2
 3L2 − 0 ⋅ L3 
w ⋅L
24⋅ E⋅ I
6 ⋅ E⋅ I
0


c =
=
=
2
24⋅ E⋅ I
4
 L3 L2 
−L


 2

 3L 2L 
So we now have the quartic equation for the deflection of the beam to be
w ⋅L
y=
0
2
2
⋅x −
w ⋅L
0
3
⋅x +
w
0
4
⋅x
24⋅ E⋅ I
12⋅ E⋅ I
24E⋅ I
If we let L = 1 and assume w = 24⋅ E⋅ I to simplify the curve, we see this graph:
0
2
1.5
2
3
4
x − 2x + x
1
0.5
0.5
0
0.5
1
1.5
x
Our interest is in the interval between 0 and 1. Recall that the graph is showing positive y direction
as up, but in our example we assumed it was down. The beam’s deflection is really as this
diagram shows:
Note that this is a little bit of a drastic deflection, but conditions here are illustrative, not
necessarily reality based. Besides, I didn't give a unit of measure, so suppose L was 100 yards.
Eigenvalues and Eigenfunctions
Suppose we desire to solve y'' + c⋅ y = 0 with boundary values y( 0 ) = 0 and y( L) = 0
Case 1 c = 0 so y'' = 0 is the DE and y = A⋅ x + B , y' = A , y'' = 0. Now y( 0 ) = 0 forces B = 0 and
y( L) = 0 means A⋅ L + 0 = 0 so A = 0 and we only have the trivial solution y = 0.
2
Case 2 c < 0 From the characteristic equation m + c = 0 with c < 0 we get solution
y = A⋅ e
− c⋅ t
+ B⋅ e
−
− c⋅ t
and from the boundary conditions y( 0 ) = 0 we get A + B = 0 or
B = −A and from y( L) = 0 we get 0 = A⋅ e
− c⋅ L
+ B⋅ e
−
− c⋅ L
or substituting
0 = A⋅ e
− A⋅ e
= 2A⋅ sinh( −c⋅ L) Since L > 0 and c < 0, 2 ⋅ A⋅ sinh( −c⋅ L) = 0 force
A = 0 = −B so again we have the trivial solution y = 0
− c⋅ L
−
− c⋅ L
Case 3 c > 0 m + c = 0 yielding y = A⋅ sin( c⋅ x) + B⋅ cos( c⋅ x). Now y( 0 ) = 0 forces B = 0 and
2
y( L) = 0 means that A⋅ sin( c⋅ L) = 0. Either A = 0 and we have the trivial solution again, OR
2
 n ⋅ π  for n = 1 , 2 , 3 , ..

 L 
c⋅ L = n ⋅ π for integer values of n. Solving for c gives c = 
This makes y = A⋅ sin( c⋅ x) a solution for any A if c =
π
2
2
,
4⋅ π
2
2
⋅
9⋅ π
2
2
, .. these values of c are
L
L
L
called characteristic values or eigenvalues for the equation y'' + c⋅ y = 0 and the corresponding
 n ⋅ π  are called characteristic functions or eigenfunctions for y'' + c⋅ y = 0
functions y = A⋅ sin

 L 
ex 1 Consider a vertical column of length L that is hinged at both ends and with a constant
vertical compressive force (pressure) P applied to its top. The boundary value problem
governing its deflection y( x) is E⋅ I⋅
d
2
y + P⋅ y = 0 with y( 0 ) = 0 and y( L) = 0 where E is a
2
dx
constant involving elasticity of the object’s material and I is the moment of inertia of a cross
section about a vertical line through its centroid. Here is a picture of the situation:
P
, y'' + c⋅ y = 0 with y( 0 ) = 0 and y( L) = 0 so
E⋅ I
 n ⋅ π ⋅ x  corresponding to eigenvalues
we either have deflection curve y = 0 or y = A⋅ sin

 L 
This is the previously described situation if c =
2
P
c =
n
n
E⋅ I
=
n ⋅π
L
2
2
for n = 1 , 2 , 3 , ..
Therefore the column will buckle only if the pressure on the column is one of the values
2 2
P =
n ⋅ π ⋅ E⋅ I
n
L
2
for n = 1 , 2 , 3 , .. and these forces are called 'critical loads'.
2
The deflection curve corresponding to the smallest critical load, called the Euler load P =
π ⋅ E⋅ I
L
2
 π ⋅ x  and is known as the first buckling mode.

 L 
is y = A⋅ sin
If a restraint is placed on the column at x =
be P =
4 ⋅ π ⋅ E⋅ I
L
2
L
2
(half-way point), then the smallest critical load will
and the deflection curve will be as in the figure shown below:
If restraints are placed at x =
L
3
and x =
2L
3
2
, the critical load becomes P =
9 ⋅ π ⋅ E⋅ I
L
and
2
so forth for other values of n.
π =0

2
ex 2 Find the eigenvalues and eigenfunctions for the BVP y'' + λ ⋅ y = 0 with y( 0 ) = 0 y'
The only non-trivial solution would be if λ > 0 (see case 1, 2, and 3 above) and then we have
the solution y = A⋅ cos( λ ⋅ x) + B⋅ sin( λ ⋅ x). The condition y( 0 ) = 0 forces A = 0 so
y' = B⋅ λ ⋅ cos( λ ⋅ x) and then and since B⋅ λ is a constant, call it C so we have
π
π
2
π
y' = C⋅ cos( λ ⋅ x) and y'  = 0 that if C ≠ 0 then λ ⋅ = ( 2n − 1 ) ⋅ so λ = ( 2n − 1 ) for
2
2
2
 
n = 1 , 2 , 3 , .. so the set of eigenvalues is λ = 1 , 9 , 25 , .. with associated eigenfunctions
y = sin( x) , y = sin( 3x) , y = sin( 5x) , ..
ex 3. The temperature u ( r) in the circular ring shown below is determined from the boundary
d
value problem: r⋅
2
2
u +
dr
u ⋅ ln
d
u = 0 with u ( a) = u and u ( b ) = u for u and u constants.
0
1
0
1
dr
r
 r
 − u 1⋅ ln a 
 
a
ln 
b
b
0
Show that u ( r) =
2 d
Multiplying both sides by r makes a Cauchy-Euler equation r ⋅
m− 1
m
solution of form u = r , u' = m⋅ r
2
(
2
)
m− 2
(
)
(
2
)
d
u + r⋅ u = 0 so we assume a
dr
dr
2
m− 2
and u'' = m − m ⋅ r
m− 1
2
so the characteristic equation is
(
m
)
2
r ⋅ m − m ⋅r
+ r⋅ m⋅ r
= 0 so we have as solution to r m − m + m = 0 a double root of
hence
the
solution
m=0
u = c + c ⋅ ln( r) (see case 2 unit III lesson 6)
1
2
Boundary condition u ( a) = u yields u = c + c ⋅ ln( a) and similarly u = c + c ⋅ ln( b ) and we can
0
solve the system:
0
1
2
1
1
2
u = c + c ⋅ ln( a)
0
1
2
u = c + c ⋅ ln( b )
1
1
2
using Cramer's rule:
 u ln( a) 
 0

 u ln( b ) 
u ⋅ ln( b ) − u ⋅ ln( a)
1
 1
 = 0
c =
1
ln( b ) − ln( a)
 1 ln( a) 


 1 ln( b) 
c =
2
Now the solution to the BVP is u = c + c ⋅ ln( r) =
1
u ⋅ ln( a) − u ⋅ ln( b )
minor adjustment u =
1
0
ln( a) − ln( b )
 1 u0 


1 u 
1

0
u −u
0
 1 ln( a) 


 1 ln( b) 
u ⋅ ln( b ) − u ⋅ ln( a)
2
+
1
ln( b ) − ln( a)
1
( ln( a) − ln( b ) )
1
1
0
ln( a) − ln( b )
0
1
0
ln( b ) − ln( a)
+
u −u
1
0
ln( b ) − ln( a)
⋅ ln( r) and a
⋅ ln( r) and adding
u ⋅ ln
u ⋅ ln( a) − u ⋅ ln( r) + u ⋅ ln( r) − u ⋅ ln( b )
u=
u −u
=
0
simplifying to u ( r) =
r
 r
 − u 1⋅ ln 
b
 a  QED
a
ln 
b