Let’s consider deflection in a simply supported (static) vertical beam: \(y(x)\), with boundary conditions \(y(0) = 0\) and \(y(L) = 0\). To get the governing equation, start with considering the sum of moments around the upper pin:
-We also know that \(M_z = E I y''\), so we can obtain
-This equation governs the stability of a beam, considering small deflections. To simplify things, let’s define \(\lambda^2 = \frac{P}{EI}\), which gives us the ODE
-We can get the general solution to this:
-To find the coefficients, let’s apply the boundary conditions, starting with \(x=0\):
-Now what? \(B \neq 0\), because otherwise we would have the trivial solution \(y(x) = 0\). Instead, to satisfy the boundary condition, we need
-\(\lambda\) give the the eigenvalues for this problem; as you can see, there are an infinite number, that correspond to eigenfunctions:
-The eigenvalues and associated eigenfunctions physically represent different modes of deflection. @@ -471,15 +471,15 @@
This means that when the combination of load force and beam properties match certain values, the beam will deflect—and buckle—in one of the modes corresponding to the associated eigenfunction.
In particular, the first mode (\(n=1\)) is interesting, because this is the first one that will be encountered if a load starts at zero and increases. This is the Euler critical load of buckling, \(P_{cr}\):
-Let’s consider a slightly different case, where at \(x=0\) the beam is supported such that \(y'(0) = 0\). How does the beam buckle in this case?
The governing equation and general solution are the same:
-but our boundary conditions are now different:
-Then, the critical buckling load, again corresponding to \(n=1\), is
-We can only get the eigenvalues analytically if we can obtain an analytical solution of the ODE, but we might want to get eigenvalues for more complex problems too. In that case, we can use an approach based on finite differences to find the eigenvalues.
Consider the same problem as above, for deflection of a simply supported beam:
-with boundary conditions \(y(0) = 0\) and \(y(L) = 0\). Let’s represent this using finite differences, for a case where \(L=3\) and \(\Delta x = 1\), so we have four points in our solution grid.
The finite difference representation of the ODE is:
-However, in this case, we are not solving for the values of deflection (\(y_i\)), but instead the eigenvalues \(\lambda\).
Then, we can write the system of equations using the above recursion formula and our two boundary conditions:
-which we can simplify down to two equations by incorporating the boundary conditions into the equations for the two middle points, and also letting \(k = \lambda^2 \Delta x^2\):
-Let’s modify this once more by multiplying both equations by \(-1\):
-Now we can represent this system of equations as a matrix equation \(A \mathbf{y} = \mathbf{b} = \mathbf{0}\):
-\(\mathbf{y} = \mathbf{0}\) is a trivial solution to this, so instead \(\det(A) = 0\) satisfies this equation. For our \(2\times 2\) matrix, that looks like:
-so the roots of this equation are \(k_1 = 1\) and \(k_2 = 3\). Recall that \(k\) is directly related to the eigenvalue: \(k = \lambda^2 \Delta x^2\), and \(\Delta x = 1\) for this case, so we can calculate the two associated eigenvalues:
-Our work has given us approximations for the first two eigenvalues. We can compare these against the exact values, given in general by \(\lambda = n \pi / L\) (which we determined above):
-So, our approximations are close, but with some obvious error. This is because we used a fairly crude step size of \(\Delta x = 1\), dividing the domain into just three segments. By using a finer resolution, we can get more-accurate eigenvalues and also more of them (remember, there are actually an infinite number!).
To do that, we will need to use Matlab, which offers the eig()
function for calculating eigenvalues—essentially it is finding the roots to the polynomial given by \(\det(A) = 0\). We need to modify this slightly, though, to use the function:
where the new matrix is
-Then, eig(A*)
will provide the values of \(k\), which we can use to find the \(\lambda\)s:
Let’s analyze the motion of masses connected by springs in a system:
- -First, we need to write the equations of motion, based on doing a free-body diagram on each mass: -\begin{align} -m_1 \frac{d^2 x_1}{dt^2} &= -k x_1 + k(x_2 - x_1) \\ -m_2 \frac{d^2 x_2}{dt^2} &= -k (x_2 - x_1) - k x_2 -\end{align} -We can condense these equations a bit: -\begin{align} -x_1^{\prime\prime} - \frac{k}{m_1} \left( -2 x_1 + x_2 \right) &= 0 \\ -x_2^{\prime\prime} - \frac{k}{m_2} \left( x_1 - 2 x_2 \right) &= 0 -\end{align} + +First, we need to write the equations of motion, based on doing a free-body diagram on each mass:
+We can condense these equations a bit:
+To proceed, we can assume that the masses will move in a sinusoidal fashion, with a shared frequency but separate amplitude:
-We can plug these into the ODEs:
-or
-Let’s put some numbers in, and try to solve for the eigenvalues: \(\omega^2\). Let \(m_1 = m_2 = 40 \) kg and \(k = 200\) N/m.
Now, the equations become
-or \(A \mathbf{y} = \mathbf{0}\), which we can represent as
-Here, \(\omega^2\) are the eigenvalues, and we can find them with \(\det(A) = 0\):
-We find there are two modes of oscillation, each associated with a different natural frequency. Unfortunately, we cannot calculate independent and unique values for the amplitudes, but if we insert the values of \(\omega\) into the above equations, we can find relations connecting the amplitudes:
-Another method of solving boundary-value problems (and also partial differential equations, as we’ll see later) involves finite differences, which are numerical approximations to exact derivatives.
Recall that the exact derivative of a function \(f(x)\) at some point \(x\) is defined as:
-So, we can approximate this derivative using a finite difference (rather than an infinitesimal difference as in the exact derivative):
-which involves some error. This is a forward difference for approximating the first derivative. We can also approximate the first derivative using a backward difference:
-To understand the error involved in these differences, we can use Taylor’s theorem to obtain Taylor series expansions:
-We can obtain higher-order approximations for the first derivative, and an approximations for the second derivative, by combining these Taylor series expansions:
-Subtracting the Taylor series for \(f(x+\Delta x)\) by that for \(f(x-\Delta x)\) gives:
-which is a second-order accurate approximation for the first derivative.
Adding the Taylor series for \(f(x+\Delta x)\) to that for \(f(x-\Delta x)\) gives:
-We can use finite differences to solve ODEs by substituting them for exact derivatives, and then applying the equation at discrete locations in the domain. This gives us a system of simultaneous equations to solve.
For example, let’s consider the ODE
-with the boundary conditions \(y(0) = 1\) and \(y(2) = 8\).
@@ -498,22 +498,22 @@To do this, we’ll follow a few steps:
1.) Replace exact derivatives in the original ODE with finite differences, and apply the equation at a particular location \((x_i, y_i)\).
For our example, this gives:
-which applies at location \((x_i, y_i)\).
2.) Next, rearrange the equation into a recursion formula:
-We can use this equation to get an equation for each of the interior points in the domain.
For the first and last points—the boundary points—we already have equations, given by the boundary conditions.
3.) Set up system of linear equations
Applying the recursion formula to the interior points, and the boundary conditions for the boundary points, we can get a system of simultaneous linear equations:
-This is a system of five equations and five unknowns, which we can solve! But, solving using substitution would be painful, so let’s represent this system of equations using a matrix and vectors:
-where \(y_1\) is the first point in the grid of points, corresponding to \(x_1 = 0\). (We saw this in the example above.) In Matlab, we can implement this equation with
@@ -647,8 +647,8 @@When a boundary condition involves a derivative, we can use a central difference to approximate the first derivative; this is more accurate than a forward or backward difference.
Consider the formula for a central difference at \(x=0\), applied for the boundary condition \(y^{\prime}(0) = 0\):
-where \(y_0\) is an imaginary, or ghost, node outside the domain. We can’t actually keep this point in our implementation, because it isn’t a real point.
We still need an equation for the point at the boundary, \(y_1\). To get this, we’ll apply the regular recursion formula, normally used at interior points:
-where \(a\), \(b\), \(c\), and \(f(x)\) depend on the problem. Normally we wouldn’t use this at the boundary node, \(y_1\), because it references a point outside the domain to the left—but we have an equation for that! From above, based on the boundary condition, we have \(y_0 = y_2\). If we incorporate that into the recursion formula, we can eliminate the ghost node \(y_0\):
-So far we’ve seen how to handle a linear boundary value problem, but what if we have a nonlinear BVP? This is going to be trickier, because our work so far relies on using linear algebra to solve the system of (linear) equations.
For example, consider the 2nd-order ODE
-with the boundary conditions \(y(0) = y(1) = 0\). This is nonlinear, due to the \(y^3\) term on the right-hand side.
To solve this, let’s first convert it into a discrete form, by replacing the second derivative with a finite difference and any \(x\)/\(y\) present with \(x_i\) and \(y_i\). We’ll also move any constants (i.e., terms that don’t contain \(y_i\)) and the nonlinear term to the right-hand side:
-where the boundary conditions are now \(y_1 = 0\) and \(y_n = 0\), with \(n\) as the number of grid points. We can rearrange and simplify into our recursion formula:
-The question is: how do we solve this now? The nonlinear term involving \(y_i^3\) on the right-hand side complicates things, but we know how to set up and solve this without the nonlinear term. We can use an approach known as successive iteration:
@@ -859,37 +859,33 @@Let’s now consider a more complicated example: heat transfer through an extended surface (a fin).
- +In this situation, we have the temperature of the body \(T_b\), the temperature of the ambient fluid \(T_{\infty}\); the length \(L\), width \(w\), and thickness \(t\) of the fin; the thermal conductivity of the fin material \(k\); and convection heat transfer coefficient \(h\).
The boundary conditions can be defined in different ways, but generally we can say that the temperature of the fin at the wall is the same as the body temperature, and that the fin is insulated at the tip. This gives us
-Our goal is to solve for the temperature distribution \(T(x)\). To do this, we need to set up a governing differential equation. Let’s do a control volume analysis of heat transfer through the fin:
- +Given a particular volumetric slice of the fin, we can define the heat transfer rates of conduction through the fin and convection from the fin to the air:
-where \(P\) is the perimeter (so that \(P \, dx\) is the heat transfer area to the fluid) and \(A_c\) is the cross-sectional area.
Performing a balance through the control volume:
-then we have as a governing equation
-where \(m^2 = (h P)/(k A_c)\).
We can obtain an exact solution for this ODE. For convenience, let’s define a new variable, \(\theta\), which is a normalized temperature:
-where \(\theta^{\prime} = T^{\prime}\) and \(\theta^{\prime\prime} = T^{\prime\prime}\). This gives us a new governing equation:
-This is a 2nd-order homogeneous ODE, which looks a lot like \(y^{\prime\prime} + a y = 0\). The exact solution is then
-We’ll use this to look at the accuracy of a numerical solution, but we will not be able to find an exact solution for more complicated versions of this problem.
We can also solve this numerically using the finite difference method. Let’s replace the derivative with a finite difference:
-which we can rearrange into a recursion formula:
-This gives us an equation for all the interior nodes; we can use the above boundary conditions to get equations for the boundary nodes. For the boundary condition at \(x=L\), \(T^{\prime}(x=L) = 0\), let’s use a backward difference:
-Let’s now consider a more-complicated case, where we also have radiation heat transfer occuring along the length of the fin. Now, our governing ODE is
-This is a bit trickier to solve because of the nonlinear term involving \(T^4\). But, we can handle it via the iterative solution method discussed above.
diff --git a/docs/content/bvps/shooting-method.html b/docs/content/bvps/shooting-method.html index 2c0cda8..9a9cf54 100644 --- a/docs/content/bvps/shooting-method.html +++ b/docs/content/bvps/shooting-method.html @@ -392,8 +392,8 @@Boundary-value problems are also ordinary differential equations—the difference is that our two constraints are at boundaries of the domain, rather than both being at the starting point.
For example, consider the ODE
-with the boundary conditions \(y(0)=1\) and \(y(2)=8\).
@@ -406,8 +406,8 @@Now, this algorithm will not work particularly well if all your guesses are random/uninformed. Fortunately, we can use linear interpolation to inform a third guess based on two initial attempts:
-Let’s try solving the given ODE using the shooting method:
-with the boundary conditions \(y(0)=1\) and \(y(2)=8\).
First, we need to convert this 2nd-order ODE into a system of two 1st-order ODEs, where we can define \(u = y'\):
-We can use the shooting method to solve a famous fluids problem: the Blasius boundary layer.
- +To get to a solveable ODE, we start with the conservation of momentum equation (i.e., Navier–Stokes equation) in the \(x\)-direction:
-and the conservation of mass equation:
-where \(u\) is the velocity component in the \(x\)-direction, \(v\) is the velocity component in the \(y\)-direction, and \(\nu\) is the fluid’s kinematic viscosity. The boundary conditions are that \(u = v = 0\) at \(y=0\), and that \(u = U_{\infty}\) as \(y \rightarrow \infty\), where \(U_{\infty}\) is the free-stream velocity.
Blasius solved this problem by converting the PDE into an ODE, by recognizing that the boundary layer thickness is given by \(\delta(x) \sim \sqrt{\frac{x \nu}{U_{\infty}}}\), and then nondimensionalizing the position coordinates using a similarity variable
-By introducing the stream function, \(\psi (x,y)\), we can ensure the continuity equation is satisfied:
-Let’s check this, using SymPy:
@@ -546,13 +544,13 @@which relates directly to the velocity components:
-with the boundary conditions \(f = f^{\prime} = 0\) at \(\eta = 0\), and \(f^{\prime} = 1\) as \(\eta \rightarrow \infty\).
This is a 3rd-order ODE, which we can solve by converting it into three 1st-order ODEs:
-with the initial condition \(y(0) = 1\). As we will see, this ODE can cause explicit numerical schemes to become unstable, and thus it is a stiff ODE. (Note that we can easily obtain the exact solution for this problem, which is \(y(t) = e^{-3 t}\).)
@@ -472,8 +472,8 @@Compare this behavior to that for the ODE
-which is non-stiff:
@@ -541,13 +541,13 @@The Backward Euler method is very similar to the Forward Euler method, except in one way: it uses the slope at the next time step:
-Then, the resulting recursion formula is
-For example, consider the problem
-To actually solve this problem with the Backward Euler method, we need to incorporate the derivative function \(f(x,y)\) into the recursion formula and solve for \(y_{i+1}\):
-This matches nearly what we saw with the Forward Euler method before—Backward Euler is also a first-order method, so the global error should be proportional to \(\Delta x\).
Let’s now return to the stiff ODE \(y^{\prime} = -3 y\), and see how the Backward Euler method does. First, we need to obtain our useable recursion formula:
-We can perform a stability analysis of the stiff problem to identify when the Forward Euler method becomes unstable. Let’s apply the method to the ODE at hand:
-where \(\sigma\) is the amplification factor. This defines whether the solution grows or decays each step—for a stable physical system, we expect the solution to get smaller or remain contant with each step.
Therefore, for the method to remain stable, we must have \(\sigma | \leq 1\). We can use this stability criterion to find conditions on \(\Delta t\) for stability:
-We can also perform a stability analysis on the Backward Euler method to show that its stability does not depend on the step size:
-For stability, we need \(| \sigma | \leq 1\):
-The temperature distribution \(T(r)\) in an annular fin of inner radius \(r_1\) and outer radius \(r_2\) is described by the equation
-where \(r\) is the radial distance from the centerline (the independent variable) and \(m^2\) is a constant that depends on the heat transfer coefficient, thermal conductivity, and thickness of the annulus. Assuming we choose a spatial step size \(\Delta r\),
- +a.) Write the finite-difference representation of the ODE (that applies at a location \(r_i\)), using central differences.
b.) Based on the last part, write the recursion formula.
c.) The boundary condition at the outer radius \(r = r_2\) is described by convection heat transfer:
-Write the boundary condition at \(r = r_2\) in recursion form (i.e., the equation you would implement into your system of equations to solve for temperature).
a.) Replace the derivatives in the given ODE with finite differences, and replace any locations with the \(i\) location:
-or
-b.) Rearrange and combine terms:
-c.) We can use a backward difference to approximate the \(dT/dr\) term. \(T_n\) represents the temperature at node \(n\) where \(r_n = r_2\):
-a.)
-The eigenfunction is then the solution function associated with an eigenvalue:
-b.) The principal eigenvalue is just that associated with \(n = 1\):
-Use the shooting method to solve the boundary value problem
-where \(y(0) = 0\) and \(y(1) = 3\). Find the initial value of \(y'\) (meaning, \(y'(0)\)) that satisfies the given boundary conditions. Use the forward Euler method with a step size of \(\Delta x = 0.5\).
First decompose into two 1st-order ODEs:
-with BCs \(z_1 (x=0) = z_{1,1} = 0\) and \(z_1(x=1) = z_{1,3} = 3\), we do not know \(y'(0) = z_2(x=0) = z_{2,1} = ?\)
Try some guess #1: \(y' (0) = 0 = z_2 (0)\), with the forward Euler method:
-so for solution 1: \(y(1) = 0 \neq 3\).
For guess #2: \(y' (0) = 2 = z_2 (0)\), with the forward Euler method:
-so for solution 1: \(y(1) = 2 \neq 3\).
For guess #3, we can interpolate:
-then, use this guess:
-so for solution 3: \(y(1) = 3\) which is the target.
So our answer is \(y'(0) = 3\).
-
-