| |
- A bar of length one is supported by an elastic underground with
spring constant c and supports a load p, nonnegative of course. Its material properties are
given by the stiffness coefficient, i.e. elasticity modulus times area inertia IE, strictly positive.
Its deformation y is given by the principle of minimizing the sum of
potential and kinetic energy, this is the sum
∫01( (1/2)(IE(x)y''(x)2+c(x)y(x)2)+p(x)y(x))dx .
By means of variational analysis one arrives at a boundary problem of a fourth order differential equation.
With no further conditions on the support one gets the natural boundary conditions y''(0)=y'''(0)=y''(1)=y'''(1)=0
given below under position 3. This requires c > 0 on an interval of nonzero length in [0,1].
We also consider the case of fixing the bar at its ends, which allows for c=0.
The variable x denotes the center line of the bar.
The differential equation reads
(-IE(x)y''(x))'' - c(x)y(x) = p(x), 0 <= x <= 1 .
- 4 correctly formulated boundary conditions make this problem uniquely solvable.
The cases implemented here are the following
-
y(0) = y'(0) = 0 and y(1) =y'(1) = 0 ,
(the bar is fixed at both ends, for example in a wall)
- y(0) = y'(0) = 0 and y''(1) = y'''(1) = 0
( the bar is fixed with zero slope on the left end and on the right
it is supported free of moments.)
- y''(0) = y'''(0) = 0 and y''(1) = y'''(1) = 0
(the bar is supported free of moments at both ends. These are the
natural boundary conditions for this case.)
In this case c must be strictly positive (on a set of nonzero measure, strictly speaking), otherwise
the problem is not well defined.
- We solve this problem by a method of finite differences.
We do not differentiate the above explicitly to obtain a fourth order
differential equation, but rather use two steps:
with
v(x) = IE(x)y''(x)
we first use the formula of order 2:
v''(xi) = (1/h2)(v(xi+1)-2v(xi)+v(xi-1)) - (1/12)h2v(4)(xi)+O(h4)
and then we substitute the same formula for y''(xj), j=i-1,i,i+1 and get
(IE(x)y''(x))''x=xj = (1/h4)(
IE(xj-1)yj-2-2(IE(xj-1)+IE(xj))yj-1
+(IE(xj-1)+4IE(xj)+IE(xj+1))yj
-2(IE(xj)+IE(xj+1))yj+1+IE(xj+1)yj+2)+O(h2).
The matrix part obtained this way is symmetric.
In order to implement the boundary conditions in a simple way we use an
equidistant grid which is shifted to the left by half the grid size.
Depending on the boundary conditions we use on the left and the right
one or two fictitious points. Hence
we use grid points xi=(i-1/2)*h, i=(-1),0,...,n+1,(n+2) with
h=1/n. There are n interior points with the numbers 1 ,..,n.
We approximate the boundary conditions as
- y1 + y0 =0, y0 - y1 =0
leading to y0 = y1 = 0 and
yn+1 + yn =0, yn+1 - yn =0
leading to yn+1 = yn = 0.
- y1 + y0 =0, y0 - y1 =0,
(difference and midvalue at x=0 are zero, both y1,y0 are zero)
and
yn+2-2yn+1+yn=0 , yn+1-2yn+yn-1=0 ,
that means that the second order approximations for y'' at 1-h/2 and
1+h/2 are both zero and hence their midvalue and their difference i.e.
y'' and y''' at x=1 too.
- y1-2y0+y-1 = 0 , y2-2y1+y0 = 0
and
yn+2-2yn+1+yn=0 , yn+1-2yn+yn-1=0 .
- This results in a linear system of equations with a banded structure with 5 diagonals.
A typical row for the interior nodes with numbers 2 to n-1 and constant IE is
-(1/h4)IE ( 0,....,0,1,-4,6,-4,1,0,.....0) - (0,....,0,c(xi),0,...,0).
The boundary conditions add two blocks two rows each on the top and the bottom of this structure.
For the conditions ''free of support'' we use the discretization of the ode also at the nodes 0 and
n+1 in order to get the necessary additional equations for the additional unknowns.
The matrix has a condition number of O(n4).
- The linear system is then solved by a special solver for banded matrices using the lu decomposition with
column pivoting.
The inverse of the matrix is componentwise positive and bounded independend of n.
We make no use of the possible simplifications by eliminating some of the variables at the boundaries
but solve the system always as one in n+2, n+3 , n+4 unknowns.
- The method is of order 2, that means that the error is of the form
e(x)*h2 (plus higher order terms)
with a smooth function e(x). By varying n and thereby h
you can study this. This however assumes that the solution y is six times continuously
differentiable.
- For larger n rounding errors in the linear systems solver may show up disastrous effects and
it makes little sense to decrease the stepsize below the fourth root of the machine precision.
Hence other solution methods using two second order equations instead are often advocated for.
|