<- previous    index    next ->

Lecture 25, Ordinary Differential Equations

An "ordinary" differential equation has exactly one independent
variable and at least one derivative of that variable.

For notation we use y=f(x) to  say y is a function of x.
For the derivative of y with respect to x, we may write
d f(x)/dx  or we can use  dy/dx  or  y'  when the independent
variable x is understood.

One of the simplest ordinary differential equations is

             y' = y

which has the analytic solution y(x) = exp(x)  the exponential function
also written as e^x .  

For a second derivative we use notation d^2 y/dx^2 or y''.

Many ordinary differential equations have analytic solutions. Yet,
many do not. We are interested in ordinary differential equations
that do not have analytic solutions and must be approximated by
numerical solutions. For testing purposes, we will start with an
equation that has an analytic solution so that we can check our
numeric solution method. The test case is the classic beam problem.

given a beam of length L, from 0 < x < L
attached rigidly at both ends
with Young's Modulus of E
with moment of inertia I
with p(x) being the load density  e.g. force per unit length
with both ends fixed, meaning:
      y=0 at x=0, y=0 at x=L, dy/dx=0 at x=0, dy/dx=0 at x=L
then the differential equation that defines the y position of
the beam at every x coordinate is

      d^4 y
   EI ----- = p(x)   with the convention for downward is negative
      dx^4

   for uniformly distributed force (weight) p(x) becomes  - W/L
   This simple case can be integrated and solved analytically:

      d^4 y 
   EI ----- = -W/L
      dx^4

     d^3 y
  EI ----- = -W/L x + A   ( A is constant of integration, value found later)
     dx^3

     d^2 y
  EI ----- = -W/L x^2/2 + A x + B 
     dx^2

     dy
  EI -- = -W/L x^3/6 + A x^2/2 + B x + C  we see C=0 because dy/dx=0 at x=0
     dx

  EI y = -W/L x^4/24 + A x^3/6 + B x^2/2 + D  we see D=0 because y=0 at x=0

  substituting y=0 at x=L in the equation above gives

     0 = -W/L L^4/24 + A L^3/6 + B L^2/2

  substituting dy/dx=0 at x=L in the equation above the above gives

     0 = -W/L L^3/6 + A L^2/2 + B L

  solving two equations in two unknowns  A = W/2  B = - WL/12
  then substituting for A and B in  EI y = ... and combining terms

         W
  y = --------  x^2 (x-L)^2
      24 L EI

  The known solution for a specific case is valuable to check your
  programming of a numerical solution before computing a more
  general case of p(x) where no closed form solution may exists.

An aside: Much information on physics equations and units and
units conversion may be found in 
units.shtml, physics equations near end

In order to find a numerical solution, set up a system of linear
equations such that the solution of the equations are values
of y(x) at some points. The codes below just use seven points, n=7,
yet could use a much larger number. Because we have a fourth order
differential equation, we will use fourth order difference equations.
The solution is fourth order (or less, actually second order) thus
we will get results with no truncation error and only small
roundoff error. The seven points will be equally spaced by h.
h = L/(n+1). The value of y(0)=0 is given and the value y(L)=0
is given. The unknown points are y(h), y(2h), ... , y(7h).

We need the difference equations for fourth order derivatives from
nderiv.out

From nderiv.out the difference equation for d y^4/dx^4
using 5 y values is:

d y^4/dx^4(x) =(1/h^4)( y(x)-4y(x+h)+6y(x+2h)-4y(x+3h)+y(x+4h) )

We will start x at h rather than zero and solve at 7 points.
Since d y^4/dx^4 = -W/L for uniform load, we can get equations
at x=2h ... x=6h. We will need special equations at x=h
and x=7h that take into account the boundary conditions..

  the partial system of equations we want looks like:

  |                             |   | y( h) |   |  0   |
  | -4   6  -4   1   0   0   0  |   | y(2h) |   | -W/L |
  |  1  -4   6  -4   1   0   0  |   | y(3h) |   | -W/L |
  |  0   1  -4   6  -4   1   0  | * | y(4h) | = | -W/L |
  |  0   0   1  -4   6  -4   1  |   | y(5h) |   | -W/L |
  |  0   0   0   1  -4   6  -4  |   | y(6h) |   | -W/L |
  |                             |   | y(7h) |   |  0   |

The boundary conditions dy/dx(0)=0, dy/dx(L)=0
must now be applied. Looking again at nderiv.out for a
first derivative that uses 5 points we find:

dy/dx(x) = 1/12h ( -25y(x) +48y(x+h) -36y(x+2h) +16y(x+3h) -3y(x+4h) )

For x=0, and knowing y(0)=0 we get the first row of the matrix

  | 48 -36  16  -3   0   0   0  |   | y( h) |   |  0   |

At x=L we find

dy/dx(x) = 1/12h ( 3y(x-4h) -16y(x-3h) +36y(x-2h) -48y(x-h) +25y(x) )

For x=L, working backward,  and knowing y(L)=0 we get
the last row of the system of equations.

  | 0   0   0  -3  16 -36  48  |   | y(7h) |   |  0   |

Solve the equations and we have y(h), y(2h), ..., y(7h) as
shown in the output files below in beam2_c.out.

In the code h=1 so that y(h) = y(x=1), y(7h) = y(x=7)
Note that the rather large step size works because we are
using a high enough order method. The analytic solution
is shown with the slopes at each point. The numeric
solution at the end checks very close to the analytic
solution.

beam2.c
beam2_c.out
beam2.f90
beam2_f90.out
beam2.adb
beam2_ada.out
beam2.java
beam2_java.out
beam2.py
beam2_py.out


For more versions of static and dynamic beam equations
Lecture 28d near the end


   <- previous    index    next ->

Other links

Go to top