This chapter focuses on methods for solving 2nd-order ODEs constrained by boundary conditions: boundary-value problems (BVPs).
- -"Eigenvalue" means characteristic value. These types of problems show up in many areas involving boundary-value problems, where we may not be able to obtain an analytical solution, but we can identify certain characteristic values that tell us important information about the system: the eigenvalues.
-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: -\begin{align} -\sum M &= M_z + P y = 0 \\ -M_z &= -P y -\end{align}
-We also know that $M_z = E I y''$, so we can obtain -\begin{align} -M_z = E I \frac{d^2 y}{dx^2} &= -P y \\ -y'' + \frac{P}{EI} y &= 0 -\end{align} -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 -\begin{equation} -y'' + \lambda^2 y = 0 -\end{equation} -We can get the general solution to this: -\begin{equation} -y(x) = A \cos (\lambda x) + B \sin (\lambda x) -\end{equation}
-To find the coefficients, let's apply the boundary conditions, starting with $x=0$: -\begin{align} -y(x=0) &= 0 = A \cos 0 + B \sin 0 \\ -\rightarrow A &= 0 \\ -y(x=L) &= 0 = B \sin (\lambda L) -\end{align} -Now what? $B \neq 0$, because otherwise we would have the trivial solution $y(x) = 0$. Instead, to satisfy the boundary condition, we need -\begin{align} -B \neq 0 \rightarrow \sin (\lambda L) &= 0 \\ -\text{so} \quad \lambda L &= n \pi \quad n = 1, 2, 3, \ldots, \infty \\ -\lambda &= \frac{n \pi}{L} \quad n = 1, 2, 3, \ldots, \infty -\end{align} -$\lambda$ give the the eigenvalues for this problem; as you can see, there are an infinite number, that correspond to eigenfunctions: -\begin{equation} -y_n = B \sin \left( \frac{n \pi x}{L} \right) \quad n = 1, 2, 3, \ldots, \infty -\end{equation}
-The eigenvalues and associated eigenfunctions physically represent different modes of deflection. -For example, consider the first three modes (corresponding to $n = 1, 2, 3$):
- -clear all; clc
-
-L = 1.0;
-x = linspace(0, L);
-subplot(1,3,1);
-y = sin(pi * x / L);
-plot(y, x); title('n = 1')
-subplot(1,3,2);
-y = sin(2 * pi * x / L);
-plot(y, x); title('n = 2')
-subplot(1,3,3);
-y = sin(3* pi * x / L);
-plot(y, x); title('n = 3')
-
-Here we see different modes of how the beam will buckle. How do we know when this happens?
-Recall that the eigenvalue is connected to the physical properties of the beam: -\begin{gather} -\lambda^2 = \frac{P}{EI} \rightarrow \lambda = \sqrt{\frac{P}{EI} = \frac{n \pi}{L} \\ -P = \frac{EI}{L} n^2 \pi^2 -\end{gather} -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}$: -\begin{gather} -\lambda_1 = \frac{\pi}{L} \rightarrow \lambda_1^2 = \frac{P}{EI} = \frac{\pi^2}{L^2} \\ -P_{cr} = \frac{\pi^2 E I}{L^2} -\end{gather}
- -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: -\begin{align} -y'' + \lambda^2 y &= 0 \\ -y(x) &= A \cos (\lambda x) + B \sin (\lambda x) -\end{align} -but our boundary conditions are now different: -\begin{align} -y'(0) = 0 = -\lambda A \sin(0) + \lambda B\cos(0) \\ -\rightarrow B &= 0 \\ -y &= A \cos (\lambda x) \\ -y(L) &= 0 = A \cos (\lambda L) \\ -A \neq 0 \rightarrow \cos(\lambda L) &= 0 \\ -\text{so} \quad \lambda L &= \frac{(2n-1) \pi}{2} \quad n = 1,2,3,\ldots, \infty \\ -\lambda &= \frac{(2n-1) \pi}{2 L} \quad n = 1,2,3,\ldots, \infty -\end{align}
-Then, the critical buckling load, again corresponding to $n=1$, is -\begin{equation} -P_{cr} = \frac{\pi^2 EI}{4 L^2} -\end{equation}
- -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: -\begin{equation} -y'' + \lambda^2 y = 0 -\end{equation} -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: -\begin{align} -\frac{y_{i-1} - 2y_i + y_{i+1}}{\Delta x^2} + \lambda^2 y_i &= 0 \\ -y_{i-1} + \left( \lambda^2 \Delta x^2 - 2 \right) y_i + y_{i+1} &= 0 -\end{align} -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: -\begin{align} -y_1 &= 0 \\ -y_1 + y_2 \left( \lambda^2 \Delta x^2 - 2 \right) + y_3 &= 0 \\ -y_2 + y_3 \left( \lambda^2 \Delta x^2 - 2 \right) + y_4 &= 0 \\ -y_4 &= 0 -\end{align} -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$: -\begin{align} -y_2 (k-2) + y_3 &= 0 \\ -y_2 + y_3 (k-2) &= 0 -\end{align} -Let's modify this once more by multiplying both equations by $-1$: -\begin{align} -y_2 (2-k) - y_3 &= 0 \\ --y_2 + y_3 (2-k) &= 0 -\end{align}
-Now we can represent this system of equations as a matrix equation $A \mathbf{y} = \mathbf{b} = \mathbf{0}$: -\begin{equation} -\begin{bmatrix} 2-k & -1 \\ -1 & 2-k \end{bmatrix} -\begin{bmatrix} y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} -\end{equation} -$\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: -\begin{align} -\det(A) = \begin{vmatrix} 2-k & -1 \\ -1 & 2-k \end{vmatrix} = (2-k)^2 - 1 &= 0 \\ -k^2 - 4k + 3 &= 0 \\ -(k-3)(k-1) &= 0 -\end{align} -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: -\begin{align} -k_1 &= \lambda_1^2 \Delta x^2 = 1 \rightarrow \lambda_1 = 1 \\ -k_2 &= \lambda_2^2 \Delta x^2 = 3 \rightarrow \lambda_2 = \sqrt{3} = 1.732 -\end{align}
-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): -\begin{align} -n=1: \quad \lambda_1 &= \frac{\pi}{L} = \frac{\pi}{3} = 1.0472 \\ -n=2: \quad \lambda_2 &= \frac{2\pi}{L} = \frac{2\pi}{3} = 2.0944 -\end{align} -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:
-\begin{align}
-\det(A) &= 0 \\
-\det \left( A^* - k I \right) = 0
-\end{align}
-where the new matrix is
-\begin{equation}
-A^* = \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix}
-\end{equation}
-Then, eig(A*) will provide the values of $k$, which we can use to find the $\lambda$s:
clear all; clc
-
-dx = 1.0;
-L = 3.0;
-
-Astar = [2 -1; -1 2];
-k = eig(Astar);
-
-lambda = sqrt(k) / dx;
-
-fprintf('lambda_1: %6.3f\n', lambda(1));
-fprintf('lambda_2: %6.3f', lambda(2));
-As expected, this matches with our manual calculation above. But, we might want to calculate these eigenvalues more accurately, so let's generalize this a bit and then try using $\Delta x= 0.1$:
- -clear all; clc
-
-dx = 0.1;
-L = 3.0;
-x = 0 : dx : L;
-n = length(x) - 2;
-
-Astar = zeros(n,n);
-for i = 1 : n
- if i == 1
- Astar(1,1) = 2;
- Astar(1,2) = -1;
- elseif i == n
- Astar(n,n-1) = -1;
- Astar(n,n) = 2;
- else
- Astar(i,i-1) = -1;
- Astar(i,i) = 2;
- Astar(i,i+1) = -1;
- end
-end
-k = eig(Astar);
-
-lambda = sqrt(k) / dx;
-
-fprintf('lambda_1: %6.3f\n', lambda(1));
-fprintf('lambda_2: %6.3f\n\n', lambda(2));
-
-err = abs(lambda(1) - (pi/L)) / (pi/L);
-fprintf('Error in lambda_1: %5.2f%%\n', 100*err);
-Let's analyze the motion of masses connected by springs in a system:
-
-
To proceed, we can assume that the masses will move in a sinusoidal fashion, with a shared frequency but separate amplitude: -\begin{align} -x_i &= A_i \sin (\omega t) \\ -x_i^{\prime\prime} &= -A_i \omega^2 \sin (\omega t) -\end{align} -We can plug these into the ODEs: -\begin{align} -\sin (\omega t) \left[ \left( \frac{2k}{m_1} - \omega^2 \right) A_1 - \frac{k}{m_1} A_2 \right] &= 0 \\ -\sin (\omega t) \left[ -\frac{k}{m_2} A_1 + \left( \frac{2k}{m_2} - \omega^2 \right) A_2 \right] &= 0 -\end{align} -or -\begin{align} -\left( \frac{2k}{m_1} - \omega^2 \right) A_1 - \frac{k}{m_1} A_2 &= 0 \\ --\frac{k}{m_2} A_1 + \left( \frac{2k}{m_2} - \omega^2 \right) A_2 &= 0 -\end{align} -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 -\begin{align} -\left( 10 - \omega^2 \right) A_1 - 5 A_2 &= 0 \\ --5 A_1 + \left( 10 - \omega^2 \right) A_2 &= 0 -\end{align} -or $A \mathbf{y} = \mathbf{0}$, which we can represent as -\begin{equation} -\begin{bmatrix} 10-\omega^2 & -5 \\ -5 & 10-\omega^2 \end{bmatrix} -\begin{bmatrix} A_1 \\ A_2 \end{bmatrix} = -\begin{bmatrix} 0 \\ 0 \end{bmatrix} -\end{equation} -Here, $\omega^2$ are the eigenvalues, and we can find them with $\det(A) = 0$: -\begin{align} -\det(B) &= 0 \\ -\det (B^* - \omega^2 I) &= 0 -\end{align}
- -clear all; clc
-
-Bstar = [10 -5; -5 10];
-omega_squared = eig(Bstar);
-omega = sqrt(omega_squared);
-
-fprintf('omega_1 = %5.2f rad/s\n', omega(1));
-fprintf('omega_2 = %5.2f rad/s\n', omega(2));
-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: -\begin{align} -\omega_1: \quad A_1 &= A_2 \\ -\omega_2: \quad A_1 &= -A_2 -\end{align}
-So, for the first mode, we have the two masses moving in sync with the same amplitude. In the second mode, they move with opposite (but equal) amplitude. With the two different frequencies, they also have two different periods:
- -t = linspace(0, 3);
-subplot(1,5,1)
-plot(sin(omega(1)*t), t); hold on
-plot(0,0, 's');
-set (gca, 'ydir', 'reverse' )
-box off; set(gca,'Visible','off')
-
-subplot(1,5,2)
-plot(sin(omega(1)*t), t); hold on
-plot(0,0, 's');
-set (gca, 'ydir', 'reverse' )
-text(-2.5,-0.2, 'First mode')
-box off; set(gca,'Visible','off')
-
-subplot(1,5,4)
-plot(-sin(omega(2)*t), t); hold on
-plot(0,0, 's');
-set (gca, 'ydir', 'reverse' )
-box off; set(gca,'Visible','off')
-
-subplot(1,5,5)
-plot(sin(omega(2)*t), t); hold on
-plot(0,0, 's');
-set (gca, 'ydir', 'reverse' )
-box off; set(gca,'Visible','off')
-text(-2.7,-0.2, 'Second mode')
-
-We can confirm that the system would actually behave in this way by setting up the system of ODEs and integrating based on initial conditions matching the amplitudes of the two modes.
-For example, let's use $x_1 (t=0) = x_2(t=0) = 1$ for the first mode, and $x_1(t=0) = 1$ and $x_2(t=0) = -1$ for the second mode. We'll use zero initial velocity for both cases.
-Then, we can solve by converting the system of two 2nd-order ODEs into a system of four 1st-order ODEs:
- -%%file masses.m
-function dxdt = masses(t, x)
-% this is a function file to calculate the derivatives associated with the system
-
-m1 = 40;
-m2 = 40;
-k = 200;
-
-dxdt = zeros(4,1);
-
-dxdt(1) = x(2);
-dxdt(2) = (k/m1)*(-2*x(1) + x(3));
-dxdt(3) = x(4);
-dxdt(4) = (k/m2)*(x(1) - 2*x(3));
-clear all; clc
-
-% this is the integration for the system in the first mode
-[t, X] = ode45('masses', [0 3], [1.0 0.0 1.0 0.0]);
-subplot(1,5,1)
-plot(X(:,1), t);
-ylabel('displacement (m)'); xlabel('time (s)')
-set (gca, 'ydir', 'reverse' )
-%box off; set(gca,'Visible','off')
-
-subplot(1,5,2)
-plot(X(:,3), t); xlabel('time (s)')
-set (gca, 'ydir', 'reverse' )
-text(-4,-0.2, 'First mode')
-
-% this is the integration for the system in the second mode
-[t, X] = ode45('masses', [0 3], [1.0 0.0 -1.0 0.0]);
-subplot(1,5,4)
-plot(X(:,1), t);
-ylabel('displacement (m)'); xlabel('time (s)')
-set (gca, 'ydir', 'reverse' )
-%box off; set(gca,'Visible','off')
-
-subplot(1,5,5)
-plot(X(:,3), t); xlabel('time (s)')
-set (gca, 'ydir', 'reverse' )
-text(-4,-0.2, 'Second mode')
-
-This shows that we get either of the pure modes of motion with the appropriate initial conditions.
-What about if the initial conditions don't match either set of amplitude patterns?
- -[t, X] = ode45('masses', [0 3], [0.25 0.0 0.75 0.0]);
-subplot(1,5,1)
-plot(X(:,1), t);
-%plot(0,0, 's');
-set (gca, 'ydir', 'reverse' )
-%box off; set(gca,'Visible','off')
-
-subplot(1,5,2)
-plot(X(:,3), t);
-set (gca, 'ydir', 'reverse' )
-
-In this case, the resulting motion will be a complicated superposition of the two modes.
- -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: -\begin{equation} -f^{\prime}(x) = \frac{df}{dx}(x) = \lim_{\Delta x \rightarrow 0} \frac{f(x+\Delta x) - f(x)}{\Delta x} -\end{equation}
-So, we can approximate this derivative using a finite difference (rather than an infinitesimal difference as in the exact derivative): -\begin{equation} -f^{\prime}(x) \approx \frac{f(x+\Delta x) - f(x)}{\Delta x} -\end{equation} -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: -\begin{equation} -f^{\prime}(x) \approx \frac{f(x) - f(x - \Delta x)}{\Delta x} -\end{equation}
-To understand the error involved in these differences, we can use Taylor's theorem to obtain Taylor series expansions: -\begin{align} -f(x + \Delta x) &= f(x) + \Delta x \, f^{\prime}(x) + \Delta x^2 \frac{1}{2!} f^{\prime\prime}(x) + \cdots \\ -\rightarrow \frac{f(x + \Delta x) - f(x)}{\Delta x} &= f^{\prime}(x) + \mathcal{O}\left( \Delta x \right) \\ -f(x - \Delta x) &= f(x) - \Delta x \, f^{\prime}(x) + \Delta x^2 \frac{1}{2!} f^{\prime\prime}(x) + \cdots \\ -\rightarrow \frac{f(x) - f(x - \Delta x)}{\Delta x} &= f^{\prime}(x) + \mathcal{O}\left( \Delta x \right) \\ -\end{align} -where the $\mathcal{O}()$ notation stands for "order of magnitude of". So, we can see that each of these approximations is first-order accurate.
- -We can obtain higher-order approximations for the first derivative, and an approximations for the second derivative, by combining these Taylor series expansions: -\begin{align} -f(x + \Delta x) &= f(x) + \Delta x \, f^{\prime}(x) + \Delta x^2 \frac{1}{2!} f^{\prime\prime}(x) + \mathcal{O}\left( \Delta x^3 \right) \\ -f(x - \Delta x) &= f(x) - \Delta x \, f^{\prime}(x) + \Delta x^2 \frac{1}{2!} f^{\prime\prime}(x) + \mathcal{O}\left( \Delta x^3 \right) -\end{align}
-Subtracting the Taylor series for $f(x+\Delta x)$ by that for $f(x-\Delta x)$ gives: -\begin{align} -f(x + \Delta x) - f(x - \Delta x) &= 2 \Delta x \, f^{\prime}(x) + \mathcal{O}\left( \Delta x^3 \right) \\ -f^{\prime}(x) &= \frac{f(x + \Delta x) - f(x - \Delta x)}{2 \Delta x} + \mathcal{O}\left( \Delta x^2 \right) -\end{align} -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: -\begin{align} -f(x + \Delta x) + f(x - \Delta x) &= 2 f(x) + \Delta x^2 f^{\prime\prime}(x) + \mathcal{O}\left( \Delta x^3 \right) \\ -f^{\prime\prime}(x) &= \frac{f(x + \Delta x) - 2 f(x) + f(x - \Delta x)}{\Delta x^2} + \mathcal{O}\left( \Delta x^2 \right) -\end{align} -which is a second-order accurate approximation for the second derivative.
- -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 -\begin{equation} -y^{\prime\prime} + x y^{\prime} - x y = 2 x \;, -\end{equation} -with the boundary conditions $y(0) = 1$ and $y(2) = 8$.
-First, we discretize the continuous domain: divide it into a number of discrete segments. For now, let's choose $\Delta x = 0.5$, which creates four segments and thus five points: $x_1 = 0, x_2 = 0.5, x_3 = 1.0, x_4 = 1.5, x_5 = 2.0$.
-Our goal is then to find approximate values of $y(x)$ at these points: $y_1$ through $y_5$. So, we have five unknowns, and need five equations to solve for them. We can use the ODE to provide these equations, by replacing the derivatives with finite differences, and applying the equation at particular discrete locations.
-Recall that $y(x)$ is a function just like $f(x)$, and so we can apply the above finite difference equations to $y(x)$ and $y(x+\Delta x)$. Now that we have points, or nodes, at locations separated by $\Delta x$, we can consider a point $x_i$ where $y(x_i) = y_i$, $y(x_i + \Delta x) = y(x_{i+1}) = y_{i+1}$, and $y(x_i - \Delta x) = y(x_{i-1}) = y_{i-1}$.
-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: -\begin{equation} -\frac{y_{i+1} - 2y_i + y_{i-1}}{\Delta x^2} + x_i \left( \frac{y_{i+1} - y_{i-1}}{2 \Delta x}\right) - x_i y_i = 2 x_i -\end{equation} -which applies at location $(x_i, y_i)$.
-2.) Next, rearrange the equation into a recursion formula: -\begin{equation} -y_{i-1} \left(1 - x_i \frac{\Delta x}{2}\right) + y_i \left( -2 -\Delta x^2 x_i \right) + y_{i+1} \left(1 + x_i \frac{\Delta x}{2}\right) = 2 x_i \Delta x^2 -\end{equation} -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: -\begin{align} -y_1 &= 1 \\ -y_1 (0.875) + y_2 (-2.125) + y_3 (1.125) &= 0.25 \\ -y_2 (0.75) + y_3 (-2.25) + y_4 (1.25) &= 0.5 \\ -y_3 (0.625) + y_4 (-2.375) + y_5 (1.375) &= 0.75 \\ -y_5 &= 8 -\end{align}
-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: -\begin{equation} -\begin{bmatrix} -1 & 0 & 0 & 0 & 0 \\ -0.875 & -2.125 & 1.125 & 0 & 0 \\ -0 & 0.75 & -2.25 & 1.25 & 0 \\ -0 & 0 & 0.625 & -2.375 & 1.375 \\ -0 & 0 & 0 & 0 & 1 -\end{bmatrix} -\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{bmatrix} = -\begin{bmatrix} 1 \\ 0.25 \\ 0.5 \\ 0.75 \\ 8 \end{bmatrix} -\end{equation} -or, more compactly, $A \mathbf{y} = \mathbf{b}$.
-4.) Solve the linear system of equations
-The final step is just to solve. We can do this in Matlab with y = A \ b. (This is equivalent to y = inv(A)*b, but faster.)
A = [1.0 0 0 0 0;
- 0.875 -2.125 1.125 0 0;
- 0 0.75 -2.25 1.25 0;
- 0 0 0.625 -2.375 1.375;
- 0 0 0 0 1];
-
-b = [1.0; 0.25; 0.5; 0.75; 8.0];
-
-x = [0 : 0.5 : 2];
-y = A \ b;
-plot(x, y, 'o-');
-
-clear all; clc
-
-dx = 0.5;
-x = [0 : dx : 2];
-n = length(x);
-A = zeros(n,n); b = zeros(n,1);
-
-for i = 1 : n
- if i == 1
- A(1,1) = 1;
- b(1) = 1;
- elseif i == n
- A(n,n) = 1;
- b(n) = 8;
- else
- A(i, i-1) = 1 - x(i)*dx/2;
- A(i, i) = -2 - x(i)*dx^2;
- A(i, i+1) = 1 + x(i)*dx/2;
- b(i) = 2*x(i)*dx^2;
- end
-end
-y = A \ b;
-plot(x, y, 'o-')
-
-This looks good, but we can get a more-accurate solution by reducing our step size $\Delta x$:
- -clear all; clc
-
-dx = 0.001;
-x = [0 : dx : 2];
-n = length(x);
-A = zeros(n,n); b = zeros(n,1);
-
-for i = 1 : n
- if i == 1
- A(1,1) = 1;
- b(1) = 1;
- elseif i == n
- A(n,n) = 1;
- b(n) = 8;
- else
- A(i, i-1) = 1 - x(i)*dx/2;
- A(i, i) = -2 - x(i)*dx^2;
- A(i, i+1) = 1 + x(i)*dx/2;
- b(i) = 2*x(i)*dx^2;
- end
-end
-y = A \ b;
-plot(x, y)
-
-We will encounter four main kinds of boundary conditions. Consider the ODE $y^{\prime\prime} + y = 0$, on the domain $0 \leq x \leq L$.
-Whichever type of boundary condition we are dealing with, the goal will be to construct an equation representing the boundary condition to incorporate in our system of equations.
-If we have a fixed value boundary condition, such as $y(0) = a$, then this equation is straightforward: -\begin{equation} -y_1 = a -\end{equation} -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
-A(1,1) = 1;
-b(1) = a;
-If we have a fixed derivative boundary condition, such as $y^{\prime}(0) = 0$, then we need to use a finite difference to represent the derivative. When the boundary condition is at the starting location, $x=0$, the easiest way to do this is with a forward difference: -\begin{align} -y^{\prime}(0) \approx \frac{y_2 - y_1}{\Delta x} &= 0 \\ --y_1 + y_2 &= 0 -\end{align} -We can implement this in Matlab with
-A(1,1) = -1;
-A(1,2) = 1;
-b(1) = 0;
-When we have this sort of derivative boundary condition at the right side of the domain, at $x=L$, then we can use a backward difference to represent the derivative: -\begin{align} -y^{\prime}(L) \approx \frac{y_n - y_{n_1}}{\Delta x} &= 0 \\ --y_{n-1} + y_n &= 0 -\end{align} -where $y_n$ is the final point ($x_n = L$) and $y_{n-1}$ is the second-to-last point ($x_{n-1} = L - \Delta x$). We can implement this in Matlab with
-A(n,n-1) = -1;
-A(n,n) = 1;
-b(n) = 0;
-If we have a linear combination of a fixed value and fixed derivative, like $a \, y(0) + b \, y^{\prime}(0) = c$, then we can combine the above approaches using a forward difference: -\begin{align} -a y(0) + b y^{\prime}(0) \approx a y_1 + b \frac{y_2 - y_1}{\Delta x} &= c \\ -(a \Delta x - b) y_1 + b y_2 &= c \Delta x -\end{align} -and in Matlab:
-A(1,1) = a*dx - b;
-A(1,2) = b;
-b(1) = c * dx;
-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$: -\begin{align} -y^{\prime}(0) \approx \frac{y_2 - y_0}{2 \Delta x} &= 0 \\ -y_0 &= y_2 -\end{align} -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: -\begin{align} -a y_{i-1} + b y_i + c y_{i+1} = f(x_i) \\ -a y_0 + b y_1 + c y_2 = f(x_1) \;, -\end{align} -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$: -\begin{align} -a y_2 + b y_1 + c y_2 &= f(x_1) \\ -b y_1 + (a + c) y_2 &= f(x_1) \, -\end{align} -which is the equation we can actually use at the boundary point. -In Matlab, this looks like
-A(1,1) = b;
-A(1,2) = a + c;
-b(1) = f(x(1));
-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 -\begin{equation} -y^{\prime\prime} = 3y + x^2 + 100 y^2 -\end{equation} -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: -\begin{equation} -\frac{y_{i-1} - 2y_i + y_{i+1}}{\Delta x^2} - 3y_i = x_i^2 + 100 y_i^2 -\end{equation} -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: -\begin{equation} -y_{i-1} + y_i \left( -2 - 3 \Delta x^2 \right) + y_{i+1} = x_i^2 \Delta x^2 + 100 \Delta x^2 y_i^2 -\end{equation}
-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:
-Let's implement that process in Matlab:
- -%% Initial setup
-clear all; clc
-
-dx = 0.01;
-x = 0 : dx : 1;
-n = length(x);
-
-A = zeros(n, n);
-b = zeros(n, 1);
-
-%% First, solve the problem without the nonlinear term:
-for i = 1 : n
- if i == 1 % x = 0 boundary condition
- A(1,1) = 1;
- b(1) = 0;
- elseif i == n % x = L boundary condition
- A(n,n) = 1;
- b(n) = 0;
- else % interior nodes, use recursion formula
- A(i, i-1) = 1;
- A(i, i) = -2 - 3*dx^2;
- A(i, i+1) = 1;
- b(i) = x(i)^2 * dx^2;
- end
-end
-% get solution without nonlinear term
-y = A \ b;
-
-plot(x, y, '--'); hold on
-
-%% Now, set up iterative process to solve while incorporating nonlinear terms
-iter = 1;
-y_old = zeros(n, 1);
-while max(abs(y - y_old)) > 1e-6
- y_old = y;
- % A matrix is not changed, but the b vector does
- for i = 2 : n - 1
- b(i) = x(i)^2 * dx^2 + 100*(dx^2)*(y_old(i)^2);
- end
-
- y = A \ b;
- iter = iter + 1;
-end
-
-fprintf('Number of iterations: %d\n', iter);
-plot(x, y)
-legend('Initial solution', 'Final solution')
-
-Another option is just to set our "guess" for the $y$ solution to be zero, rather than solve the problem in two steps:
- -%% Initial setup
-clear all; clc
-
-dx = 0.01;
-x = 0 : dx : 1;
-n = length(x);
-
-A = zeros(n, n);
-b = zeros(n, 1);
-
-%% Set up the coefficient matrix, which does not change
-for i = 1 : n
- if i == 1 % x = 0 boundary condition
- A(1,1) = 1;
- b(1) = 0;
- elseif i == n % x = L boundary condition
- A(n,n) = 1;
- b(n) = 0;
- else % interior nodes, use recursion formula
- A(i, i-1) = 1;
- A(i, i) = -2 - 3*dx^2;
- A(i, i+1) = 1;
- b(i) = x(i)^2 * dx^2;
- end
-end
-
-% just use zeros as our initial guess for the solution
-y = zeros(n, 1);
-
-%% Successive iteration
-iter = 1;
-y_old = 100 * rand(n, 1); % setting this to some random values, just to enter the while loop
-while max(abs(y - y_old)) > 1e-6
- y_old = y;
- % A matrix is not changed, but the b vector does
- for i = 2 : n - 1
- b(i) = x(i)^2 * dx^2 + 100*(dx^2)*(y_old(i)^2);
- end
-
- y = A \ b;
- iter = iter + 1;
-end
-
-fprintf('Number of iterations: %d\n', iter);
-This made our process take slightly more iterations, because the initial guess was slightly further away from the final solution. For other problems, having a bad initial guess could make the process take much longer, so coming up with a good initial guess may be important.
- -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 -\begin{align} -T(x=0) &= T_b \\ -q(x=L) = 0 \rightarrow \frac{dT}{dx} (x=0) &= 0 -\end{align}
-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: -\begin{align} -q_{\text{conv}} &= h P \left( T - T_{\infty} \right) dx \\ -q_{\text{cond}, x} &= -k A_c \left(\frac{dT}{dx}\right)_{x} \\ -q_{\text{cond}, x+\Delta x} &= -k A_c \left(\frac{dT}{dx}\right)_{x+\Delta x} \;, -\end{align} -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: -\begin{align} -q_{\text{cond}, x+\Delta x} &= q_{\text{cond}, x} - q_{\text{conv}} \\ --k A_c \left(\frac{dT}{dx}\right)_{x+\Delta x} &= -k A_c \left(\frac{dT}{dx}\right)_{x} - h P \left( T - T_{\infty} \right) dx \\ --k A_c \frac{\left.\frac{dT}{dx}\right|_{x+\Delta x} - \left.\frac{dT}{dx}\right|_{x}}{dx} &= -h P ( T - T_{\infty} ) \\ -\lim_{\Delta x \rightarrow 0} : -k A_c \left. \frac{d^2 T}{dx^2} \right|_x &= -h P (T - T_{\infty}) \\ -\frac{d^2 T}{dx^2} &= \frac{h P}{k A_c} (T - T_{\infty}) \\ -\frac{d^2 T}{dx^2} &= m^2 (T - T_{\infty}) -\end{align} -then we have as a governing equation -\begin{equation} -\frac{d^2 T}{dx^2} - m^2 (T - T_{\infty}) = 0 \;, -\end{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: -\begin{equation} -\theta \equiv T - T_{\infty} -\end{equation} -where $\theta^{\prime} = T^{\prime}$ and $\theta^{\prime\prime} = T^{\prime\prime}$. -This gives us a new governing equation: -\begin{equation} -\theta^{\prime\prime} - m^2 \theta = 0 \;. -\end{equation} -This is a 2nd-order homogeneous ODE, which looks a lot like $y^{\prime\prime} + a y = 0$. The exact solution is then -\begin{align} -\theta(x) &= c_1 e^{-m x} + c_2 e^{m x} \\ -T(x) &= T_{\infty} + c_1 e^{-m x} + c_2 e^{m x} -\end{align} -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: -\begin{align} -\frac{d^2 T}{dx^2} - m^2 (T - T_{\infty}) &= 0 \\ -\frac{T_{i-1} - 2T_i + T_{i+1}}{\Delta x^2} - m^2 \left( T_i - T_{\infty} \right) &= 0 -\end{align} -which we can rearrange into a recursion formula: -\begin{equation} -T_{i-1} + T_i \left( -2 - \Delta x^2 m^2 \right) + T_{i+1} = -m^2 \Delta x^2 \, T_{\infty} -\end{equation} -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: -\begin{align} -T_1 &= T_b \\ -\frac{T_n - T_{n-1}}{\Delta x} = 0 \rightarrow - T_{n-1} + T_n &= 0 -\end{align}
-Combining all these equations, we can construct a linear system: $A \mathbf{T} = \mathbf{b}$.
- -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 -\begin{equation} -\frac{d^2 T}{dx^2} - \frac{h P}{k A_c} \left(T - T_{\infty}\right) - \frac{\sigma \epsilon P}{h A_c} \left(T^4 - T_{\infty}^4 \right) = 0 -\end{equation}
-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.
- -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 -\begin{equation} -y^{\prime\prime} + xy^{\prime} - xy = 2x -\end{equation} -with the boundary conditions $y(0)=1$ and $y(2)=8$.
-The numerical methods we have already discussed (e.g., Forward Euler, Runge-Kutta) require values of $y$ and $y^{\prime}$ at the starting point, $x=0$. So we can't use these directly because we are missing $y^{\prime}(0)$.
-But, what if we could guess a value for the missing initial condition, then integrate towards the second boundary condition using one of our familiar numerical methods, and then adjust our guess if necessary and repeat? This concept is the shooting method.
-The shooting method algorithm is:
-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: -\begin{align} -\text{guess 3} &= \text{guess 2} + m \left( \text{target} - \text{solution 2} \right) \\ -m &= \frac{\text{guess 1} - \text{guess 2}}{\text{solution 1} - \text{solution 2}} -\end{align} -where "target" is the target boundary condition—in this case, $y(L)$.
- -Let's try solving the given ODE using the shooting method: -\begin{equation} -y^{\prime\prime} + xy^{\prime} - xy = 2x -\end{equation} -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'$: -\begin{align} -y' &= u \\ -u' &= 2x + xy - xu -\end{align}
- -%%file shooting_rhs.m
-function dydx = shooting_rhs(x, y)
-
-dydx = zeros(2,1);
-dydx(1) = y(2);
-dydx(2) = 2*x - x*y(2) + x*y(1);
-clear all; clc
-
-% target boundary condition
-target = 8;
-
-% Pick a guess for y'(0) of 1
-guess1 = 1;
-[X, Y] = ode45('shooting_rhs', [0 2], [1 guess1]);
-solution1= Y(end,1);
-fprintf('Solution 1: %5.2f\n', solution1);
-
-% Pick a second guess for y'(0) of 4
-guess2 = 4;
-[X, Y] = ode45('shooting_rhs', [0 2], [1 guess2]);
-solution2 = Y(end,1);
-fprintf('Solution 2: %5.2f\n', solution2);
-
-% now use linear interpolation to find a new guess
-m = (guess1 - guess2)/(solution1 - solution2);
-guess3 = guess2 + m*(target-solution2);
-fprintf('Guess 3: %5.2f\n', guess3);
-
-[X, Y] = ode45('shooting_rhs', [0 2], [1 guess3]);
-solution3 = Y(end,1);
-fprintf('Solution 3: %5.2f\n', solution3);
-fprintf('Target: %5.2f\n', target);
-
-plot(X, Y(:,1)); axis([0 2 0 9])
-
-As you can see, using linear interpolation, we are able to find the correct guess for the missing initial condition $y'(0)$ with in just three steps. This works so well because this is a linear ODE. If we had a nonlinear ODE, it would take more tries, as we'll see shortly.
- -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: -\begin{equation} -u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \nu \frac{\partial^2 u}{\partial y^2} -\end{equation} -and the conservation of mass equation: -\begin{equation} -\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \;, -\end{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 -\begin{equation} -\eta = y \sqrt{\frac{U_{\infty}}{2 \nu x}} -\end{equation}
-By introducing the stream function, $\psi (x,y)$, we can ensure the continuity equation is satisfied: -\begin{equation} -u = \frac{\partial \psi}{\partial y} \;, \quad v = -\frac{\partial \psi}{\partial x} -\end{equation}
-Let's check this, using SymPy:
- -
%%python
-import sympy as sym
-sym.init_printing()
-x, y, u, v = sym.symbols('x y u v')
-
-# Streamfunction
-psi = sym.Function(r'psi')(x,y)
-
-# Define u and v based on the streamfunction
-u = psi.diff(y)
-v = -psi.diff(x)
-
-# Check the continuity equation:
-print(u.diff(x) + v.diff(y) == 0)
-Using the boundary layer thickness and free-stream velocity, we can define the dimensionlesss stream function $f(\eta)$: -\begin{equation} -f(\eta) = \frac{\psi}{U_{\infty}} \sqrt{\frac{U_{\infty}}{2 \nu x}} -\end{equation} -which relates directly to the velocity components: -\begin{align} -u &= \frac{\partial \psi}{\partial y} = \frac{\partial \psi}{\partial f} \frac{\partial f}{\partial \eta} \frac{\partial \eta}{\partial y} \\ - &= U_{\infty} \sqrt{\frac{2 \nu x}{U_{\infty}}} \cdot f^{\prime}(\eta) \cdot \sqrt{\frac{U_{\infty}}{2 \nu x}} \\ -u &= U_{\infty} f^{\prime} (\eta) \\ -v &= -\frac{\partial \psi}{\partial x} = -\left( \frac{\partial \psi}{\partial x} + \frac{\partial \psi}{\partial \eta} \frac{\partial \eta}{\partial x} \right) \\ - &= \sqrt{\frac{\nu U_{\infty}}{2x}} \left( \eta f^{\prime} - f \right) -\end{align}
-We can insert these into the $x$-momentum equation, which leads to an ODE for the dimensionless stream function $f(\eta)$: -\begin{equation} -f^{\prime\prime\prime} + f f^{\prime\prime} = 0 \;, -\end{equation} -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: -\begin{align} -y_1 &= f \quad y_1^{\prime} = y_2 \\ -y_2 &= f^{\prime} \quad y_2^{\prime} = y_3 \\ -y_3 &= f^{\prime\prime} \quad y_3^{\prime} = -y_1 y_3 -\end{align} -and we can use the shooting method to solve by recognizing that we have two initial conditions, $y_1(0) = y_2(0) = 0$, and are missing $y_3(0)$. We also have a target boundary condition: $y_2(\infty) = 1$.
-(Note: obviously we cannot truly integrate over $0 \leq \eta < \infty$. Instead, we just need to choose a large enough number. In this case, using 10 is sufficient.)
-Let's create a function to evaluate the derivatives:
- -%%file blasius_rhs.m
-function dydx = blasius_rhs(eta, y)
-
-dydx = zeros(3,1);
-
-dydx(1) = y(2);
-dydx(2) = y(3);
-dydx(3) = -y(1) * y(3);
-First, let's try the same three-step approach we used for the simpler example, taking two guesses and then using linear interpolation to find a third guess:
- -clear all; clc
-
-target = 1.0;
-
-guesses = zeros(3,1);
-solutions = zeros(3,1);
-
-guesses(1) = 1;
-[eta, F] = ode45('blasius_rhs', [0 10], [0 0 guesses(1)]);
-solutions(1) = F(end, 2);
-
-guesses(2) = 0.1;
-[eta, F] = ode45('blasius_rhs', [0 10], [0 0 guesses(2)]);
-solutions(2) = F(end, 2);
-
-m = (guesses(1) - guesses(2))/(solutions(1) - solutions(2));
-guesses(3) = guesses(2) + m*(target - solutions(2));
-
-[eta, F] = ode45('blasius_rhs', [0 10], [0 0 guesses(3)]);
-solutions(3) = F(end, 2);
-
-tries = [1; 2; 3];
-table(tries, guesses, solutions)
-fprintf('Target: %5.2f\n', target);
-So, for this problem, using linear interpolation did not get us the correct solution on the third try. This is because the ODE is nonlinear. But, you can see that we are converging towards the correct solution—it will just take more tries.
-Rather than manually take an unknown (and potentially large) number of guesses, let's automate this with a while loop:
clear all; clc
-
-target = 1.0;
-
-% get these arrays of stored values started.
-% note: I'm only doing this to make it easier to show a table of values
-% at the end; otherwise, there's no need to store these values.
-tries = [1; 2; 3];
-guesses = zeros(3,1);
-solutions = zeros(3,1);
-
-guesses(1) = 1;
-[eta, F] = ode45('blasius_rhs', [0 10], [0 0 guesses(1)]);
-solutions(1) = F(end, 2);
-
-guesses(2) = 0.1;
-[eta, F] = ode45('blasius_rhs', [0 10], [0 0 guesses(2)]);
-solutions(2) = F(end, 2);
-
-num = 2;
-solutions(3) = -1000.; % doing this to kick off the while loop
-while abs(target - solutions(num)) > 1.e-9
- num = num + 1;
- m = (guesses(num-2) - guesses(num-1))/(solutions(num-2) - solutions(num-1));
- guesses(num) = guesses(num-1) + m*(target - solutions(num-1));
- [eta, F] = ode45('blasius_rhs', [0 1e3], [0 0 guesses(num)]);
- solutions(num) = F(end, 2);
- tries(num) = num;
-
- % we should probably set a maximum number of iterations, just to prevent
- % an infinite while loop in case something goes wrong
- if num >= 1e4
- break
- end
-end
-
-table(tries, guesses, solutions)
-fprintf('Number of iterations required: %d', num)
-%plot -r 200
-plot(F(:, 2), eta); ylim([0 5])
-xlabel("f^{\prime}(\eta) = u/U_{\infty}")
-ylabel('\eta')
-
-We can see that this plot of $\eta$, the $y$ position normalized by the boundary-layer thickness, vs. nondimensional velocity matches the original figure.
- -Welcome to the jupyter-book repository! We're excited you're here and want to contribute.
These guidelines are designed to make it as easy as possible to get involved. -If you have any questions that aren't discussed below, please let us know by opening an issue!
-Before you start you'll need to set up a free GitHub account and sign in. -Here are some instructions.
-Already know what you're looking for in this guide? Use the TOC to the right -to navigate this page!
-jupyter-book is a young project maintained by a growing group of enthusiastic developers— and we're excited to have you join!
-Most of our discussions will take place on open issues.
As a reminder, we expect all contributors to jupyter-book to adhere to the Jupyter Code of Conduct in these conversations.
git is a really useful tool for version control. -GitHub sits on top of git and supports collaborative and distributed working.
-You'll use Markdown to chat in issues and pull requests on GitHub.
-You can think of Markdown as a few little symbols around your text that will allow GitHub
-to render the text with formatting.
-For example you could write words as bold (**bold**), or in italics (*italics*),
-or as a link ([link](https://https://youtu.be/dQw4w9WgXcQ)) to another webpage.
GitHub has a helpful page on -getting started with writing and formatting Markdown on GitHub.
-Every project on GitHub uses issues slightly differently.
-The following outlines how the jupyter-book developers think about these tools.
Issues are individual pieces of work that need to be completed to move the project forwards. -A general guideline: if you find yourself tempted to write a great big issue that -is difficult to describe as one unit of work, please consider splitting it into two or more issues.
-Issues are assigned labels which explain how they relate to the overall project's -goals and immediate next steps.
-The current list of labels are here and include:
-
These issues contain a task that a member of the team has determined we need additional help with.
If you feel that you can contribute to one of these issues, we especially encourage you to do so!
-
These issues contain a task that a member of the team thinks could be a good entry point to the project.
If you're new to the jupyter-book project, we think that this is a great place for your first contribution!
These issues point to problems in the project.
If you find new a bug, please give as much detail as possible in your issue, including steps to recreate the error. - If you experience the same bug as one already listed, please add any additional information that you have as a comment.
-
These issues are asking for enhancements to be added to the project.
Please try to make sure that your enhancement is distinct from any others that have already been requested or implemented. - If you find one that's similar but there are subtle differences please reference the other request in your issue.
-
These are questions that users and contributors have asked.
Please check the issues (especially closed ones) to see if your question has been asked and answered before. - If you find one that's similar but there are subtle differences please reference the other request in your issue.
-This section covers the general structure of the -Jupyter Book repository, and -explains which pieces are where.
-The Jupyter Book repository contains two main pieces:
-This is used to help create and build books.
-It can be found at ./jupyter_book.
page module builds single pages. This module is meant to be self-contained for
-converting single .ipynb/.md/etc pages into HTML. Jupyter Book uses this module when
-building entire books, but the module can also be used on its own (it's what jupyter-book page uses).
-You can find the module at: jupyter_book/page.create.py and build.py create and build a book. They connect with the CLI and
-are used to process multiple pages and stitch them together into a static website template.This is used when generating new books. This website defines the structure of
-the site that is created when you run jupyter-book create. It contains the Javascript, CSS, and
-HTML structure of a book. It can be found at
-jupyter_book/book_template.
_includes/
-folder contains core HTML and javascript files for the site. For example,
-_includes/head.html contains the HTML for the header of each page, which is where CSS and JS files are linked.assets/
-folder contains static CSS/JS files that don't depend on site configuration._sass/
-folder contains all of the book and page CSS rules. This is stitched together in a single CSS file
-at build time (SCSS is a way to split up CSS rules among multiple files). Within this folder, the
-_sass/page/ folder
-has CSS files for a single page of content, while the other folders/files contain CSS rules for
-the whole book.content/
-folder contains the content for the Jupyter Book documentation (e.g., the markdown for this page).Here are a few examples of how this code gets used to help you get started.
-jupyter-book create mybook/, the create.py module is used to generate an empty template using the template in jupyter_book/book_template/.jupyter-book build mybook/, the build.py module to loop through your page content files,
-and uses the page/ module to convert each one into HTML and places it in mybook/_build.Hopefully this explanation gets you situated and helps you understand how the pieces all fit together. -If you have any questions, feel free to open an issue asking for help!
-We appreciate all contributions to jupyter-book, but those accepted fastest will follow a workflow similar to the following:
1. Comment on an existing issue or open a new issue referencing your addition.
-This allows other members of the jupyter-book development team to confirm that you aren't overlapping with work that's currently underway and that everyone is on the same page with the goal of the work you're going to carry out.
-This blog is a nice explanation of why putting this work in up front is so useful to everyone involved.
-2. Fork the jupyter-book repository to your profile.
-This is now your own unique copy of jupyter-book. -Changes here won't effect anyone else's work, so it's a safe space to explore edits to the code!
-Make sure to keep your fork up to date with the master repository.
-3. Make the changes you've discussed.
-Try to keep the changes focused. -We've found that working on a new branch makes it easier to keep your changes targeted.
-4. Submit a pull request.
-A member of the development team will review your changes to confirm that they can be merged into the main code base. -When opening the pull request, we ask that you follow some specific conventions. -We outline these below.
-To improve understanding pull requests "at a glance", we encourage the use of several standardized tags. -When opening a pull request, please use at least one of the following prefixes:
-You can also combine the tags above, for example if you are updating both a test and -the documentation: [TST, DOC].
-Pull requests should be submitted early and often!
-If your pull request is not yet ready to be merged, please open your pull request as a draft. -More information about doing this is available in GitHub's documentation. -This tells the development team that your pull request is a "work-in-progress", -and that you plan to continue working on it.
-When your pull request is Ready for Review, you can select this option on the PR's page, -and a project maintainer will review your proposed changes.
-We welcome and recognize all contributions from documentation to testing to code development. -You can see a list of current contributors in the contributors tab.
-You're awesome.
-— Based on contributing guidelines from the STEMMRoleModels project.
- -When equations are of this form, we can directly integrate:
-\begin{align} -\frac{dy}{dx} &= y^{\prime} = f(x) \\ -\int dy &= \int f(x) dx \\ -y(x) &= \int f(x) dx + C -\end{align}For example: -\begin{align} -\frac{dy}{dx} &= x^2 \\ -y(x) &= \frac{1}{3} x^3 + C -\end{align}
-While these problems look simple, there may not be an obvious closed-form solution to all:
-\begin{align} -\frac{dy}{dx} &= e^{-x^2} \\ -y(x) &= \int e^{-x^2} dx + C -\end{align}(You may recognize this as leading to the error function, $\text{erf}$: -$\frac{1}{2} \sqrt{\pi} \text{erf}(x) + C$, -so the exact solution to the integral over the range $[0,1]$ is 0.7468.)
-If the given derivative is a separate function of $x$ and $y$, then we can solve via separation of variables: -\begin{align} -\frac{dy}{dx} &= f(x) g(y) = \frac{h(x)}{j(y)} \\ -\int \frac{1}{g(y)} dy &= \int f(x) dx -\end{align}
-For example, consider this problem: -\begin{equation} -y^{\prime} = \frac{dy}{dx} = 1 + y^2 \\ -\end{equation} -We can separate this into a problem that looks like $f(y) dy = g(x) dx$, where $dy = \frac{1}{1+y^2}$ and $g(x) = 1$. -\begin{align} -\int \frac{dy}{1 + y^2} &= \int dx \\ -\arctan y &= x + c \\ -y(x) &= \tan(x+c) -\end{align}
-Unfortunately, not every separable ODE can be integrated: -\begin{align} -\frac{dy}{dx} &= \frac{e^x / 2 + 5}{y^2 + \cos y} \\ -(y^2 + \cos y) dy &= (e^x / 2 + 5) dx -\end{align}
-Given a general linear 1st-order ODE of the form -\begin{equation} -\frac{dy}{dx} + p(x) y = q(x) -\end{equation} -we can solve by integration factor: -\begin{equation} -y(x) = e^{-\int p(x) dx} \left[ \int e^{\int p(x) dx} q(x) dx + C \right] -\end{equation}
-For example, in this equation -\begin{equation} -y^{\prime} + xy - 5 e^x = 0 -\end{equation} -after rearranging to the standard form -\begin{equation} -y^{\prime} + xy = 5 e^x -\end{equation} -we see that $p(x) = x$ and $q(x) = 5e^x$.
-Given a general nonlinear 1st-order ODE -\begin{equation} -\frac{dy}{dx} + p(x) y = q(x) y^a -\end{equation} -where $a \neq 1$ and $a$ is a constant. This is known as the Bernoulli equation.
-We can solve by transforming to a linear equation, by changing the dependent variable from $y$ to $z$: -\begin{align} -\text{let} \quad z &= y^{1-a} \\ -\frac{dz}{dx} &= (1-a) y^{-a} \frac{dy}{dx} -\end{align} -Multiply the original equation by $(1-a) y^{-a}$: -\begin{align} -(1-a) y^{-a} \frac{dy}{dx} + (1-a) y^{-a} p(x) y &= (1-a) y^{-a} q(x) y^a \\ -\frac{dz}{dx} + p(x) (1-a) z &= q(x) (1-a) \;, -\end{align} -which is now a linear first-order ODE, that looks like -\begin{equation} -\frac{dz}{dx} + p(x)^{\prime} z = q(x)^{\prime} -\end{equation} -where $p(x)^{\prime} = (1-a) p(x)$ and $q(x)^{\prime} = (1-a)q(x)$.
-We can solve this using the integrating-factor approach discussed above. Then, once we have $z(x)$, we can find $y(x)$: -\begin{align} -z &= y^{1-a} \\ -y &= z^{\frac{1}{1-a}} -\end{align}
- -This website is an interactive Jupyter Book for ME 373, Mechanical Engineering Methods, taught at Oregon State University.
- -Applying the trapezoidal rule and Simpson's rule introduces the concept of error in numerical solutions.
-In our work so far, we have come across two obvious kinds of error, that we'll come back to later:
-In any numerical solution, there are five main sources of error:
-Error in input data: this comes from measurements, and can be systematic (for example, due to uncertainty in measurement devices) or random.
-Rounding errors: loss of significant digits. This comes from the fact that computers cannot represent real numbers exactly, and instead use a floating-point representation.
-Truncation error: due to an infinite process being broken off. For example, an infinite series or sum ending after a finite number of terms, or discretization error by using a finite step size to approximate a continuous function.
-Error due to simplifications in a mathematical model: "All models are wrong, but some are useful" (George E.P. Box) All models make some idealizations, or simplifying assumptions, which introduce some error with respect to reality. For example, we may assume gases are continuous, that a spring has zero mass, or that a process is frictionless.
-Human error and machine error: there are many potential sources of error in any code. These can come from typos, human programming errors, errors in input data, or (more rarely) a pure machine error. Even textbooks, tables, and formulas may have errors.
-We can also differentiate between absolute and relative error in a quantity. If $y$ is an exact value and $\tilde{y}$ is an approximation to that value, then we have
-If $y$ is a vector, then we can define error using the maximum of the elements: -\begin{equation} -\max_i \frac{ |y_i - \tilde{y}_{i} |}{|y_i|} -\end{equation}
- -For numerically solving definite integrals ($\int_a^b f(x) dx$) we have methods like the trapezoidal rule and Simpson's rule. When we need to solve 1st-order ODEs of the form -\begin{equation} -y^{\prime} = \frac{dy}{dx} = f(x, y) -\end{equation} -for $y(x)$, we need other methods. All of them will work by starting at the initial conditions, and then using information provided by the ODE to march forward in the solution, based on an increment (i.e., step size) $\Delta x$.
-For example, let's say we want to solve -\begin{equation} -\frac{dy}{dx} = 4 x - \frac{2 y}{x} \;, \quad y(1) = 1 -\end{equation} -This problem is fairly simple, and we can find the general and particular solutions to compare our numerical results against: -\begin{align} -\text{general: } y(x) &= x^2 + \frac{x}{x^2} \\ -\text{particular: } y(x) &= x^2 -\end{align}
-Recall that the derivative, $y^{\prime}$, is the same as the slope. At the starting point, $(x,y) = (1,1)$, where $y^{\prime} = 2$, this looks like:
- -format compact
-%plot inline
-
-x = linspace(1, 3);
-y = x.^2;
-plot(x, y); hold on
-plot([1, 2], [1, 3], '--')
-legend(['Solution'], ['Slope at start'])
-hold off
-
-Remember that the slope, or derivative, is -\begin{equation} -\text{slope} = \frac{\text{rise}}{\text{run}} = \frac{\Delta y}{\Delta x} -\end{equation}
-Let's consider the initial condition—the starting point—as $(x_i, y_i)$, and the next point in our numerical solution is $(x_{i+1}, y_{i+1})$, where $i$ represents an index starting at 1 and ending at the number of steps $N$. Our step size is then $\Delta x = x_{i+1} - x_i$.
-Based on our (very simple) approximation to the first derivative based on slope, we can relate the derivative to our two points: -\begin{equation} -\left(\frac{dy}{dx}\right)_{i} = \frac{y_{i+1} - y_i}{x_{i+1} - x_i} = \frac{y_{i+1} - y_i}{\Delta x} -\end{equation} -Then, solve this for our unknown: -\begin{equation} -y_{i+1} = y_i + \left(\frac{dy}{dx}\right)_i \Delta x -\end{equation} -This is the Forward Euler method.
-Based on a given step size $\Delta x$, we'll use this formula (called a recursion formula) to march forward and obtain the full solution over given $x$ domain. That will look something like this:
- -clear x y
-x_exact = linspace(1, 3);
-y_exact = x_exact.^2;
-plot(x_exact, y_exact); hold on
-
-% our derivative function, dy/dx
-f = @(x,y) 4*x - (2*y)/x;
-
-dx = 0.5;
-x = 1 : dx : 3;
-y(1) = 1;
-for i = 1 : length(x)-1
- y(i+1) = y(i) + f(x(i), y(i))*dx;
-end
-plot(x, y, 'o--', 'MarkerFaceColor', 'r')
-
-legend(['Exact solution'], ['Numerical solution'], 'Location','northwest')
-hold off
-
-Another way to obtain the recursion formula for the Forward Euler method is to use a Taylor series expansion. -Recall that for well-behaved functions, the Taylor series expansion says -\begin{equation} -y(x + \Delta x) = y(x) + \Delta x y^{\prime}(x) + \frac{1}{2} \Delta x^2 y^{\prime\prime}(x) + \frac{1}{3!} \Delta x^3 y^{\prime\prime\prime}(x) \dots \;. -\end{equation} -This is exact for an infinite series. We can apply this formula to our (unknown) solution $y_i$ and cut off the terms of order $\Delta x^2$ and higher; the derivative $y^{\prime}$ is given by our original ODE. -This gives us the same recursion formula as above: -\begin{equation} -\therefore y_{i+1} \approx y_i + \left( \frac{dy}{dx}\right)_i \Delta x -\end{equation} -where we can now see that we are introducing some error on the order of $\Delta x^2$ at each step. This is the local truncation error. The global error is the accumulation of error over all the steps, and is on the order of $\Delta x$. Thus, the Forward Euler method is a first-order method, because its global error is on the order of the step size to the first power: error $\sim \mathcal{O}(\Delta x)$.
-Forward Euler is also an explicit method, because its recursion formula is explicity defined for $y_{i+1}$. (You'll see when that may not be the case soon.)
-In general, for an $n$th-order method: -\begin{align} -\text{local error } &\sim \mathcal{O}(\Delta x^{n+1}) \\ -\text{global error } &\sim \mathcal{O}(\Delta x^{n}) -\end{align} -(This only applies for $\Delta x < 1$; in cases where you have a $\Delta x > 1$, you should nondimensionalize the problem based on the domain size such that $0 \leq x \leq 1$.)
-Applying the Forward Euler method then requires:
-Let's do another example: -\begin{equation} -y^{\prime} = 8 e^{-x}(1+x) - 2y -\end{equation} -with the initial condition $y(0) = 1$, and the domain $0 \leq x \leq 7$. This is a linear 1st-order ODE that we can find the analytical solution for comparison: -\begin{equation} -y(x) = e^{-2x} (8 x e^x + 1) -\end{equation}
-To solve, we'll create an anonymous function for the derivative and then incorporate that into our Forward Euler code. We'll start with $\Delta x = 0.2$.
- -clear
-
-f = @(x,y) 8*exp(-x)*(1 + x) - 2*y;
-
-dx = 0.2;
-x = 0 : dx : 7;
-n = length(x);
-y(1) = 1;
-
-% Forward Euler loop
-for i = 1 : n - 1
- y(i+1) = y(i) + dx*f(x(i), y(i));
-end
-
-x_exact = linspace(0, 7);
-y_exact = exp(-2.*x_exact).*(8*x_exact.*exp(x_exact) + 1);
-plot(x_exact, y_exact); hold on
-plot(x, y, 'o--')
-legend('Exact solution', 'Forward Euler solution')
-
-Notice the visible error in that plot, which is between 0.2–0.25, or in other words $\mathcal{O}(\Delta x)$.
-How can we reduce the error? Just like with the trapezoidal rule, we have two main options:
-The downside to reducing $\Delta x$ is the increased number of steps we then have to take, which may make the solution too computationally expensive. A more-accurate method would have less error per step, which might allow us to use the same $\Delta x$ but get a better solution. Let's next consider some better methods.
-Heun's method is a predictor-corrector method; these work by predicting a solution at some intermediate location and then using that information to get a better overall answer at the next location (correcting). Heun's uses the Forward Euler method to predict the solution at $x_{i+1}$, then uses the average of the slopes at $y_i$ and the predicted $y_{i+1}$ to get a better overall answer for $y_{i+1}$.
-\begin{align} -\text{predictor: } y_{i+1}^p &= y_i + \Delta x f(x_i, y_i) \\ -\text{corrector: } y_{i+1} &= y_i + \frac{\Delta x}{2} \left( f(x_i, y_i) + f(x_{i+1}, y_{i+1}^p) \right) -\end{align}Heun's method is second-order accurate, meaning the global error is $\mathcal{O}(\Delta x^2)$ and explicit.
-Let's see this method in action:
- -clear
-
-f = @(x,y) 8*exp(-x)*(1 + x) - 2*y;
-
-dx = 0.2;
-x = 0 : dx : 7;
-n = length(x);
-y(1) = 1;
-
-% Heun's method loop
-for i = 1 : n - 1
- y_p = y(i) + dx*f(x(i), y(i));
- y(i+1) = y(i) + (dx/2)*(f(x(i), y(i)) + f(x(i+1), y_p));
-end
-
-x_exact = linspace(0, 7);
-y_exact = exp(-2.*x_exact).*(8*x_exact.*exp(x_exact) + 1);
-plot(x_exact, y_exact); hold on
-plot(x, y, 'o--')
-legend('Exact solution', "Heun's method solution")
-fprintf('Maximum error: %5.3f', abs(max(y_exact) - max(y)))
-
-Notice how the error is visibly smaller than for the Forward Euler method–the maximum error is around 0.05, which is very close to $\Delta x^2 = 0.04$.
-The midpoint method, also known as the modified Euler method, is another predictor-corrector method, that instead predicts the solution at the midpoint ($x + \Delta x/2$): -\begin{align} -y_{i + \frac{1}{2}} &= y_i + \frac{\Delta x}{2} f(x_i, y_i) \\ -y_{i+1} &= y_i + \Delta x f \left( x_{i+\frac{1}{2}} , y_{i + \frac{1}{2}} \right) -\end{align}
-Like Heun's method, the midpoint method is explicit and second-order accurate:
- -clear
-
-f = @(x,y) 8*exp(-x)*(1 + x) - 2*y;
-
-dx = 0.2;
-x = 0 : dx : 7;
-n = length(x);
-y(1) = 1;
-
-% midpoint method loop
-for i = 1 : n - 1
- y_half = y(i) + (dx/2)*f(x(i), y(i));
- y(i+1) = y(i) + dx * f(x(i) + dx/2, y_half);
-end
-
-x_exact = linspace(0, 7);
-y_exact = exp(-2.*x_exact).*(8*x_exact.*exp(x_exact) + 1);
-plot(x_exact, y_exact); hold on
-plot(x, y, 'o--')
-legend('Exact solution', "Midpoint method solution")
-
-fprintf('Maximum error: %5.3f', abs(max(y_exact) - max(y)))
-
-Runge–Kutta methods are a family of methods that use one or more stages; the methods we have discussed so far (Forward Euler, Heun's, and midpoint) actually all fall in this family. There is also a popular fourth-order method: the fourth-order Runge–Kutta method (RK4). This uses four stages to get a more accurate solution: -\begin{align} -y_{i+1} &= y_i + \frac{\Delta x}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) \\ -k_1 &= f(x_i, y_i) \\ -k_2 &= f \left( x_i + \frac{\Delta x}{2}, y_i + \frac{\Delta x}{2} k_1 \right) \\ -k_3 &= f \left( x_i + \frac{\Delta x}{2}, y_i + \frac{\Delta x}{2} k_2 \right) \\ -k_4 &= f \left( x_i + \Delta x, y_i + \Delta x \, k_3 \right) -\end{align}
-This method is explicit and fourth-order accurate: error $\sim \mathcal{O}(\Delta x^4)$:
- -clear
-
-f = @(x,y) 8*exp(-x)*(1 + x) - 2*y;
-
-dx = 0.2;
-x = 0 : dx : 7;
-n = length(x);
-y(1) = 1;
-
-% 4th-order Runge-Kutta method loop
-for i = 1 : n - 1
- k1 = f(x(i), y(i));
- k2 = f(x(i) + dx/2, y(i) + dx*k1/2);
- k3 = f(x(i) + dx/2, y(i) + dx*k2/2);
- k4 = f(x(i) + dx, y(i) + dx*k3);
- y(i+1) = y(i) + (dx/6) * (k1 + 2*k2 + 2*k3 + k4);
-end
-
-x_exact = linspace(0, 7);
-y_exact = @(x) exp(-2.*x).*(8*x.*exp(x) + 1);
-plot(x_exact, y_exact(x_exact)); hold on
-plot(x, y, 'o--')
-legend('Exact solution', "RK4 solution")
-
-fprintf('Maximum error: %6.4f', max(abs(y_exact(x) - y)))
-
-The maximum error (0.0004) is actually a bit smaller than $\Delta x^4 = 0.0016$, but approximately the same order of magnitude.
-Matlab also offers a built-in RK4 integrator: ode45. (It is actually slightly more complicated than the equations shown just now, because it automatically adjusts the step size $\Delta x$ to control error.) You can call this function with the syntax:
[X, Y] = ode45(function_name, [x_start x_end], [IC]);
-
-where function_name is the name of a function that provides the derivative (this can be a regular function given in a file, or an anonymous function); [x_start x_end] provides the domain of integration ($x_{\text{start}} \leq x \leq x_{\text{end}}$), and [IC] provides the initial condition $y(x=x_{\text{start}})$.
For example, let's use this and compare with our exact solution:
- -clear
-
-f = @(x,y) 8*exp(-x)*(1 + x) - 2*y;
-
-[X, Y] = ode45(f, [0 7], [1]);
-
-x_exact = linspace(0, 7);
-y_exact = @(x) exp(-2.*x).*(8*x.*exp(x) + 1);
-plot(x_exact, y_exact(x_exact)); hold on
-plot(X, Y, 'o--')
-legend('Exact solution', "ode45 solution")
-
-fprintf('Maximum error: %6.4f', max(abs(y_exact(X) - Y)))
-
-What about when we cannot integrate a function analytically? In other words, when there is no (obvious) closed-form solution. In these cases, we can use numerical methods to solve the problem.
-Let's use this problem: -\begin{align} -\frac{dy}{dx} &= e^{-x^2} \\ -y(x) &= \int e^{-x^2} dx + C -\end{align}
-(You may recognize this as leading to the error function, $\text{erf}$: -$\frac{1}{2} \sqrt{\pi} \text{erf}(x) + C$, -so the exact solution to the integral over the range $[0,1]$ is 0.7468.)
- -x = linspace(0, 1);
-f = @(x) exp(-x.^2);
-plot(x, f(x))
-axis([0 1 0 1])
-
-In such cases, we can find the integral by using the trapezoidal rule, which finds the area under the curve by creating trapezoids and summing their areas: -\begin{equation} -\text{area under curve} = \sum \left( \frac{f(x_{i+1}) + f(x_i)}{2} \right) \Delta x -\end{equation}
-Let's see what this looks like with four trapezoids ($\Delta x = 0.25$):
- -hold off
-x = linspace(0, 1);
-plot(x, f(x)); hold on
-axis([0 1 0 1])
-
-x = 0 : 0.25 : 1;
-
-% plot the trapezoids
-for i = 1 : length(x)-1
- xline = [x(i), x(i)];
- yline = [0, f(x(i))];
- line(xline, yline, 'Color','red','LineStyle','--')
- xline = [x(i+1), x(i+1)];
- yline = [0, f(x(i+1))];
- line(xline, yline, 'Color','red','LineStyle','--')
- xline = [x(i), x(i+1)];
- yline = [f(x(i)), f(x(i+1))];
- line(xline, yline, 'Color','red','LineStyle','--')
-end
-hold off
-
-Now, let's integrate using the trapezoid formula given above:
- -dx = 0.1;
-x = 0.0 : dx : 1.0;
-
-area = 0.0;
-for i = 1 : length(x)-1
- area = area + (dx/2)*(f(x(i)) + f(x(i+1)));
-end
-
-fprintf('Numerical integral: %f\n', area)
-exact = 0.5*sqrt(pi)*erf(1);
-fprintf('Exact integral: %f\n', exact)
-fprintf('Error: %f %%\n', 100.*abs(exact-area)/exact)
-We can see that using the trapezoidal rule, a numerical integration method, with an internal size of $\Delta x = 0.1$ leads to an approximation of the exact integral with an error of 0.08%.
-You can make the trapezoidal rule more accurate by:
-First, how does reducing the segment size (step size) by a factor of 10 affect the error?
- -dx = 0.01;
-x = 0.0 : dx : 1.0;
-
-area = 0.0;
-for i = 1 : length(x)-1
- area = area + (dx/2)*(f(x(i)) + f(x(i+1)));
-end
-
-fprintf('Numerical integral: %f\n', area)
-exact = 0.5*sqrt(pi)*erf(1);
-fprintf('Exact integral: %f\n', exact)
-fprintf('Error: %f %%\n', 100.*abs(exact-area)/exact)
-So, reducing our step size by a factor of 10 (using 100 segments instead of 10) reduced our error by a factor of 100!
- -We can increase the accuracy of our numerical integration approach by using a more sophisticated interpolation scheme with each segment. In other words, instead of using a straight line, we can use a polynomial. Simpson's rule, also known as Simpson's 1/3 rule, refers to using a quadratic polynomial to approximate the line in each segment.
-Simpson's rule defines the definite integral for our function $f(x)$ from point $a$ to point $b$ as -\begin{equation} -\int_a^b f(x) \approx \frac{1}{6} \Delta x \left( f(a) + 4 f \left(\frac{a+b}{2}\right) + f(b) \right) -\end{equation} -where $\Delta x = b - a$.
-That equation comes from interpolating between points $a$ and $b$ with a third-degree polynomial, then integrating by parts.
- -hold off
-x = linspace(0, 1);
-plot(x, f(x)); hold on
-axis([-0.1 1.1 0.2 1.1])
-
-plot([0 1], [f(0) f(1)], 'Color','black','LineStyle',':');
-
-% quadratic polynomial
-a = 0; b = 1; m = (b-a)/2;
-p = @(z) (f(a).*(z-m).*(z-b)/((a-m)*(a-b))+f(m).*(z-a).*(z-b)/((m-a)*(m-b))+f(b).*(z-a).*(z-m)/((b-a)*(b-m)));
-plot(x, p(x), 'Color','red','LineStyle','--');
-
-xp = [0 0.5 1];
-yp = [f(0) f(m) f(1)];
-plot(xp, yp, 'ok')
-hold off
-legend('exact', 'trapezoid fit', 'polynomial fit', 'points used')
-
-We can see that the polynomial fit, used by Simpson's rule, does a better job of of approximating the exact function, and as a result Simpson's rule will be more accurate than the trapezoidal rule.
-Next let's apply Simpson's rule to perform the same integration as above:
- -dx = 0.1;
-x = 0.0 : dx : 1.0;
-
-area = 0.0;
-for i = 1 : length(x)-1
- area = area + (dx/6.)*(f(x(i)) + 4*f((x(i)+x(i+1))/2.) + f(x(i+1)));
-end
-
-fprintf('Simpson rule integral: %f\n', area)
-exact = 0.5*sqrt(pi)*erf(1);
-fprintf('Exact integral: %f\n', exact)
-fprintf('Error: %f %%\n', 100.*abs(exact-area)/exact)
-Simpson's rule is about three orders of magnitude (~1000x) more accurate than the trapezoidal rule.
-In this case, using a more-accurate method allows us to significantly reduce the error while still using the same number of segments/steps.
- -This chapter describes numerical methods used to solve integrals and first-order ordinary differential equations, along with concepts related to these such as error and stability.
- -In the past when you've talked about stability, it has likely been regarding the stability of a system. Stable systems are those will well-behaved exact solutions, meaning they do not grow unbounded. -In engineering we mostly focus (or want!) stable systems, although there are some interesting unstable systems such as those involving resonance, nonlinear dynamics, or chaos—generally we want to know when that happens so we can prevent it.
-We can also define the stability of a numerical scheme, which is when the numerical solution exhibits unphysical behavior. In other words, it blows up.
-For example, let's consider the relatively simple 1st-order ODE -\begin{equation} -\frac{dy}{dt} = -3 y -\end{equation} -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}$.)
-Let's try solving this with the Forward Euler method, integrating over $0 \leq t \leq 10$, for a range of time-step size values: $\Delta t = 0.1, 0.25, 0.5, 0.75$:
- -clear
-f = @(t,y) -3*y;
-
-dt = 0.1;
-t = 0 : dt : 20;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,1);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.25;
-t = 0 : dt : 20;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,2);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.5;
-t = 0 : dt : 20;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,3);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.75;
-t = 0 : dt : 20;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,4);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-At the smaller step sizes, $\Delta t = 0.1$ and $\Delta t = 0.25$, we see that the solution is well-behaved. But, when we increase $\Delta t$ to 0.5, we see some instability that goes away with time. Then, when we increase $\Delta t$ to 0.75, the solution eventually blows up, leading to error much larger than what we should expect based on the method's order of accuracy (first) and the step size value.
-Compare this behavior to that for the ODE -\begin{equation} -\frac{dy}{dt} = e^{-t} -\end{equation} -which is non-stiff:
- -clear
-f = @(t,y) exp(-t);
-
-dt = 0.1;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,1);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.25;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,2);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.5;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,3);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.75;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) + dt*f(t(i), y(i));
-end
-subplot(4,1,4);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-In this case, we see that the solution remains well-behaved even for larger time-step sizes, and the error matches the expected order based on the method and step-size value.
-In general numerical schemes can be:
-Schemes may be stable for some problem/system and not for another, and vice versa.
-Stability is related to robustness of a method, which is generally a tradeoff between complexity and computational cost. The choice of method and solution strategy depends on what you want, and how long you can wait for it. In general, we almost always want to use the largest $\Delta t$ allowable.
-Rather than reducing $\Delta t$ to avoid stability problems, we can also use a method that is unconditionally stable, such as the Backward Euler method.
- -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: -\begin{equation} - \left(\frac{dy}{dx}\right)_{i+1} \approx \frac{y_{i+1} - y_i}{\Delta x} -\end{equation} -Then, the resulting recursion formula is -\begin{equation} -y_{i+1} = y_i + \Delta x \left(\frac{dy}{dx}\right)_{i+1}, \text{ or} \\ -y_{i+1} = y_i + \Delta x \, f(x_{i+1}, y_{i+1}) -\end{equation} -where $f(x,y) = dy/dx$.
-Notice that this recursion formula cannot be directly solved, because $y_{i+1}$ shows up on both sides. This is an implicit method, where all the other methods we have covered (Forward Euler, Heun's, Midpoint, and 4th-order Runge-Kutta) are explicit. Implicit methods require more work to actually implement.
-For example, consider the problem -\begin{equation} -\frac{dy}{dx} = f(x,y) = 8 e^{-x} (1+x) - 2y -\end{equation} -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}$: -\begin{align} -y_{i+1} &= y_i + \Delta x \, f(x_{i+1}, y_{i+1}) \\ -y_{i+1} &= y_i + \Delta x \left[ 8 e^{-x_{i+1}} (1 + x_{i+1}) - 2 y_{i+1} \right] \\ -y_{i+1} &= y_i + 8 e^{-x_{i+1}} (1 + x_{i+1}) \Delta x - 2 \Delta x \, y_{i+1} \\ -y_{i+1} + 2 \Delta x \, y_{i+1} &= y_i + 8 e^{-x_{i+1}} (1 + x_{i+1}) \Delta x \\ -y_{i+1} &= \frac{ y_i + 8 e^{-x_{i+1}} (1 + x_{i+1}) \Delta x }{ 1 + 2 \Delta x } -\end{align} -Now we have a useable recursion formula that we can use to solve this problem. Let's use the initial condition $y(0) = 1$, the domain $0 \leq x \leq 7$, and $\Delta x = 0.2$.
- -clear
-
-dx = 0.2;
-x = 0 : dx : 7;
-y = zeros(length(x), 1);
-y(1) = 1;
-
-% Backward Euler loop
-for i = 1 : length(x) - 1
- y(i+1) = (y(i) + 8*exp(-x(i+1))*(1 + x(i+1))*dx) / (1 + 2*dx);
-end
-
-x_exact = linspace(0, 7);
-y_exact = exp(-2.*x_exact).*(8*x_exact.*exp(x_exact) + 1);
-plot(x_exact, y_exact); hold on
-plot(x, y, 'o--')
-legend('Exact solution', 'Backward Euler solution')
-
-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: -\begin{align} -y_{i+1} &= y_i + \Delta t \, f(t_{i+1}, y_{i+1}) \\ -y_{i+1} &= y_i + \Delta t \, \left( -3 y_{i+1} \right) \\ -y_{i+1} + 3 y_{i+1} \Delta t &= y_i \\ -y_{i+1} &= \frac{y_i}{1 + 3 \Delta t} -\end{align}
- -clear
-
-dt = 0.1;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) / (1 + 3*dt);
-end
-subplot(4,1,1);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.25;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) / (1 + 3*dt);
-end
-subplot(4,1,2);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.5;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) / (1 + 3*dt);
-end
-subplot(4,1,3);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-dt = 0.75;
-t = 0 : dt : 10;
-y = zeros(length(t), 1);
-y(1) = 1;
-for i = 1 : length(t) - 1
- y(i+1) = y(i) / (1 + 3*dt);
-end
-subplot(4,1,4);
-plot(t, y); title(sprintf('dt = %4.2f', dt));
-
-In this case, we see that the solution remains well-behaved for all the step sizes, not showing any of the instability we saw with the Forward Euler method. This is because the Backward Euler method is unconditionally stable.
- -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: -\begin{align} -\frac{dy}{dt} &= -3 y \\ -y_{i+1} &= y_i + \Delta t f(t_i, y_i) \\ -y_{i+1} &= y_i + \Delta t (-3 y_i) \\ - &= y_i (1 - 3 \Delta t) \\ -\frac{y_{i+1}}{y_i} &= \sigma = 1 - 3 \Delta t -\end{align} -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: -\begin{gather} -| \sigma | = | 1 - 3 \Delta t | \leq 1 \\ --1 \leq 1 - 3 \Delta t \leq 1 \\ --1 \leq 1 - 3 \Delta t \quad \text{or} \quad 1 - 3 \Delta t \leq 1 \\ -\frac{-2}{3} \leq -\Delta t \quad \quad -\Delta t \leq 0 \\ -\rightarrow \Delta t \leq \frac{2}{3} \quad \text{and} \quad \Delta t \geq 0 \\ -\therefore 0 \leq \Delta t \leq \frac{2}{3} -\end{gather} -for stability. (For safety, we might use $\Delta t < 1/2$ for safety, to stay away from the absolute stability limit.)
-The Forward Euler method is then conditionally stable.
-As a general rule of thumb, all explicit methods are conditionally stable; these are methods where the recursion formula for $y_{i+1}$ can be written and calculated explicitly in terms of known quantities.
-We can also perform a stability analysis on the Backward Euler method to show that its stability does not depend on the step size: -\begin{align} -\frac{dy}{dt} &= -3 y \\ -y_{i+1} &= y_i + \Delta t f(t_{i+1}, y_{i+1}) \\ -y_{i+1} &= y_i + \Delta t (-3 y_{i+1}) \\ -y_{i+1} (1 + 3 \Delta t) &= y_i \\ -\sigma &= \frac{y_{i+1}}{y_i} = \frac{1}{1 + 3 \Delta t} -\end{align} -For stability, we need $| \sigma | \leq 1$: -\begin{align} -| \sigma | &= \left| \frac{1}{1 + 3 \Delta t} \right| \leq 1 \\ -\rightarrow \Delta t &> 0 -\end{align} -Therefore the Backward Euler method is unconditionally stable.
- -The classic example of an elliptic PDE is Laplace's equation (yep, the same Laplace that gave us the Laplace transform), which in two dimensions for a variable $u(x,y)$ is -\begin{equation} -\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \nabla^2 u = 0 \;, -\end{equation} -where $\nabla$ is del, or nabla, and represents the gradient operator: $\nabla = \frac{\partial}{\partial x} + \frac{\partial}{\partial y}$.
-Laplace's equation shows up in a number of physical problems, including heat transfer, fluid dynamics, and electrostatics. For example, the heat equation for conduction in two dimensions is -\begin{equation} -\frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \;, -\end{equation} -where $u(x,y,t)$ is temperature and $\alpha$ is thermal diffusivity. Steady-state heat transfer (meaning after any initial transient period) is then described by Laplace's equation.
-A related elliptic PDE is Poisson's equation: -\begin{equation} -\nabla^2 u = f(x,y) \;, -\end{equation} -which also appears in multiple physical problems—most notably, when solving for pressure in the Navier–Stokes equations.
-To numerically solve these equations, and any elliptic PDE, we can use finite differences, where we replace the continuous $x,y$ domain with a discrete grid of points. This is similar to what we did with boundary-value problems in one dimension—but now we have two dimensions.
-To approximate the second derivatives in Laplace's equation, we can use central differences in both the $x$ and $y$ directions, applied around the $u_{i,j}$ point: -\begin{align} -\frac{\partial^2 u}{\partial x^2} &\approx \frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{\Delta x^2} \\ -\frac{\partial^2 u}{\partial y^2} &\approx \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{\Delta y^2} -\end{align} -where $i$ is the index used in the $x$ direction, $j$ is the index in the $y$ direction, and $\Delta x$ and $\Delta y$ are the step sizes in the $x$ and $y$ directions. -In other words, $x_i = (i-1) \Delta x$ and $y_j = (j-1) \Delta y$.
-The following figure shows the points necessary to approximate the partial derivatives in the PDE at a location $(x_i, y_j)$, for a general 2D region. This is known as a five-point stencil:
-
- Applying these finite differences gives us an approximation for Laplace's equation: -\begin{equation} -\frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{\Delta x^2} + \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{\Delta y^2} = 0 \;. -\end{equation} -If we use a uniform grid where $\Delta x = \Delta y = h$, then we can simplify to -\begin{equation} -u_{i+1,j} + u_{i,j+1} + u_{i-1,j} + u_{i,j-1} - 4 u_{i,j} = 0 \;. -\end{equation}
- -As an example, let's consider the problem of steady-state heat transfer in a square solid object. If $u(x,y)$ is temperature, then this is described by Laplace's equation: -\begin{equation} -\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \nabla^2 u = 0 \;, -\end{equation} -and we can solve this using finite differences. Using a uniform grid where $\Delta x = \Delta y = h$, Laplace's equation gives us a recursion formula that relates the values at neighboring points: -\begin{equation} -u_{i+1,j} + u_{i,j+1} + u_{i-1,j} + u_{i,j-1} - 4 u_{i,j} = 0 \;. -\end{equation}
-Consider a case where the square has sides of length $L$, and the boundary conditions are that the temperature is fixed at 100 on the left, right, and bottom sides, and fixed at 0 on the top. -For now, we'll use two segments to discretize the domain in each directions, giving us nine total points in the grid. -The following figures show the example problem, and the grid of points we'll use.
-
- 
- |
-
- 
- |
-
Using the above recursion formula, we can write an equation for each of the nine unknown points (in the interior, not the boundary points): -\begin{align} -u_{1,1} &= 100 \\ -u_{2,1} &= 100 \\ -u_{3,1} &= 100 \\ -u_{1,2} &= 100 \\ -\text{for } u_{2,2}: \quad u_{3,2} + u_{2,3} + u_{1,2} + u_{2,1} - 4u_{2,2} &= 0 \\ -u_{3,2} &= 100 \\ -u_{1,3} &= 100 \\ -u_{2,3} &= 0 \\ -u_{3,3} &= 100 -\end{align} -where $u_{i,j}$ are the unknowns. Note that in this we used the side boundary condition values for the corner points $u_{1,3}$ and $u_{3,3}$, rather than the top value. (In reality this would represent a discontinuity in temperature, so these aren't very realistic boundary conditions.)
-This is a system of linear equations, that we can represent as a matrix-vector product: -\begin{align} -\begin{bmatrix} -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ -0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\ -0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ -0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ -0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ -0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\\end{bmatrix} -\begin{bmatrix} u_{1,1} \\ u_{2,1} \\ u_{3,1} \\ u_{1,2} \\ u_{2,2} \\ u_{3,2} \\ u_{1,3} \\ u_{2,3} \\ u_{3,3} \end{bmatrix} &= -\begin{bmatrix} 100 \\ 100 \\ 100 \\ -100 \\ 0 \\ 100 \\ -100 \\ 0 \\ 100 \end{bmatrix} \\ -\text{or} \quad A \mathbf{u} &= \mathbf{b} -\end{align} -where $A$ is a $9\times 9$ coefficient matrix, $\mathbf{u}$ is a nine-element vector of unknown variables, and $\mathbf{b}$ is a nine-element right-hand side vector. -For $\mathbf{u}$, we had to take variables that physically represent points in a two-dimensional space and combine them in some order to form a one-dimensional column vector. Here, we used a row-major mapping, where we started with the point in the first row and first column, then added the remaining points in that row, before moving to the next row and repeating. We'll discuss this a bit more later.
-If we set this up in Matlab, we can solve with u = A \ b:
clear all; clc
-
-A = [
-1 0 0 0 0 0 0 0 0;
-0 1 0 0 0 0 0 0 0;
-0 0 1 0 0 0 0 0 0;
-0 0 0 1 0 0 0 0 0;
-0 1 0 1 -4 1 0 1 0;
-0 0 0 0 0 1 0 0 0;
-0 0 0 0 0 0 1 0 0;
-0 0 0 0 0 0 0 1 0;
-0 0 0 0 0 0 0 0 1];
-b = [100; 100; 100; 100; 0; 100; 100; 0; 100];
-
-% solve system of linear equations
-u = A \ b;
-
-disp(u)
-This gives us the values for temperature at each of the nine points. In this example, we really only have one unknown temperature: $u_{2,2}$, located in the middle. Does the value given make sense? We can check by rearranging the recursion formula for Laplace's equation: -\begin{equation} -u_{i,j} = \frac{u_{i+1,j} + u_{i,j+1} + u_{i-1,j} + u_{i,j-1}}{4} \;, -\end{equation} -which shows that in such problems the value of the middle point should be the average of the four surrounding points. This matches the value of 75 found above.
-We can use a contour plot to visualize the results, though we'll need to convert the one-dimensional solution array into a two-dimensional matrix to plot. The Matlab reshape() function can help us here: it reshapes an array into a matrix, by specifying the target number of desired columns and rows:
% Example of using the reshape function, with a simple array going from 1 to 10
-
-% We want to convert it into a matrix with 5 columns and 2 rows.
-% The expected output is:
-% [1 2 3 4 5;
-% 6 7 8 9 10]
-
-b = (1 : 10)';
-A = reshape(b, [5, 2]);
-disp('b array:')
-disp(b)
-disp('A matrix:')
-disp(A)
-This behavior may be a bit unexpected, because reshape() uses a column-major mapping. We can fix this by taking the transpose of the resulting matrix:
disp('transpose of output matrix:')
-disp(A')
-% We can use the reshape function to convert the calculated temperatures
-% into a 3x3 matrix:
-
-n = 3; m = 3;
-u_square = reshape(u, [n, m]);
-
-contourf(u_square')
-colorbar
-
-Overall that looks correct: the boundary conditions are right, and we see that the center is the average of the boundaries.
-But, clearly only using nine points (with eight of those being boundary conditions) doesn't give us a very good solution. To make this more accurate, we'll need to use more points, which also means we need to automate the construction of the system of equations.
- -For a two-dimensional elliptic PDE like Laplace's equation, we can generate a general recursion formula, but we need a way to take a grid of points where location is defined by row and column index and map these into a one-dimensional column vector, which has its own index.
-The following figure shows a general 2D grid of points, with $n$ number of columns in the $x$ direction (using index $i$) and $m$ number of rows in the $y$ direction (using index $j$):
-
- We want to convert the rows and columns of $u_{i,j}$ points defined by column and row index into a single column array using a different index, $k$ (this choice is arbitrary): -\begin{equation} -\begin{bmatrix} u_{1,1} \\ u_{2,1} \\ u_{3,1} \\ \vdots \\ u_{n,1} \\ -u_{1,2} \\ u_{2,2} \\ u_{3,2} \\ \vdots \\ u_{n, 2} \\ u_{1,3} \\ \vdots \\ -u_{1,m} \\ u_{2,m} \\ \vdots \\ u_{n,m} -\end{bmatrix} -\end{equation} -where $k$ refers to the index used in that array.
-To do this mapping, we can use this formula: -\begin{equation} -k_{i,j} = (j-1)n + i -\end{equation} -where $k_{i,j}$ refers to the 1D index $k$ mapped from the 2D indices $i$ and $j$.
-
- For example, in this $3\times 3$ grid, where $n=3$ and $m=3$, consider the point where $i=2$ and $j=2$ (the point right in the center). Using our formula, -\begin{equation} -k_{2,2} = (2-1)3 + 2 = 5 -\end{equation} -which matches what we can visually confirm.
-Using that mapping, we can also identify the 1D indices associated with the points surrounding location $(i,j)$: -\begin{align} -k_{i-1,j} &= (j-1)n + i - 1 \\ -k_{i+1,j} &= (j-1)n + i + 1 \\ -k_{i,j-1} &= (j-2)n + i \\ -k_{i,j+1} &= j n + i -\end{align} -which we can use to determine the appropriate locations to place values in the coefficient and right-hand side matrices.
- -clear; clc; close all
-
-h = 0.1;
-x = [0 : h : 1]; n = length(x);
-y = [0 : h : 1]; m = length(y);
-
-% The coefficient matrix A is now m*n by m*n, since that is the total number of points.
-% The right-hand side vector b is m*n by 1.
-A = zeros(m*n, m*n);
-b = zeros(m*n, 1);
-
-u_left = 100;
-u_right = 100;
-u_bottom = 100;
-u_top = 0;
-
-for j = 1 : m
- for i = 1 : n
- % for convenience we calculate all the indices once
- kij = (j-1)*n + i;
- kim1j = (j-1)*n + i - 1;
- kip1j = (j-1)*n + i + 1;
- kijm1 = (j-2)*n + i;
- kijp1 = j*n + i;
-
- if i == 1
- % this is the left boundary
- A(kij, kij) = 1;
- b(kij) = u_left;
- elseif i == n
- % right boundary
- A(kij, kij) = 1;
- b(kij) = u_right;
- elseif j == 1
- % bottom boundary
- A(kij, kij) = 1;
- b(kij) = u_bottom;
- elseif j == m
- % top boundary
- A(kij, kij) = 1;
- b(kij) = u_top;
- else
- % these are the coefficients for the interior points,
- % based on the recursion formula
- A(kij, kim1j) = 1;
- A(kij, kip1j) = 1;
- A(kij, kijm1) = 1;
- A(kij, kijp1) = 1;
- A(kij, kij) = -4;
- end
- end
-end
-u = A \ b;
-
-u_square = reshape(u, [n, m]);
-contourf(x, y, u_square')
-c = colorbar;
-c.Label.String = 'Temperature';
-
-So far, we have only discussed cases where we have Dirichlet boundary conditions; in other words, when we have all fixed values at the boundary. Frequently we also encounter Neumann-style boundary conditions, where we have the derivative specified at the boundary.
-We can handle this in the same way we do for one-dimensional boundary value problems: either with a forward or backward difference (both of which are first-order accurate), or with a central difference using an imaginary point/ghost node (which is second-order accurate). Let's focus on using the central difference, since it is more accurate.
-
- For example, let's say that at the upper boundary, the derivative of temperature is zero: -\begin{equation} -\left. \frac{\partial u}{\partial y} \right|_{\text{boundary}} = 0 -\end{equation}
-Let's consider this boundary condition applied at the point shown, $u_{2,3}$. -We can approximate this derivative using a central difference: -\begin{align} -\frac{u_{2,3}}{\partial y} \approx \frac{u_{2,4} - u_{2,2}}{\Delta x} &= 0 \\ -u_{2,4} &= u_{2,2} -\end{align} -This tells us the value of the point above the boundary, $u_{2,4}$; however, this point is a "ghost" or imaginary point located outside the boundary, so we don't really care about its value. Instead, we can use this relationship to give us a usable equation for the boundary point, by incorporating it into the normal recursion formula for Laplace's equation: -\begin{align} -u_{1,3} + u_{3,3} + u_{2,4} + u_{2,2} - 4u_{2,3} &= 0 \\ -u_{1,3} + u_{3,3} + u_{2,2} + u_{2,2} - 4u_{2,3} &= 0 \\ -\rightarrow u_{1,3} + u_{3,3} + 2 u_{2,2} - 4u_{2,3} &= 0 -\end{align}
-The recursion formula for points along the upper boundary would then become -\begin{equation} -u_{i+1,j} + u_{i-1,j} + 2 u_{i,j-1} - 4 u_{i,j} = 0 \;. -\end{equation}
-Now let's try solving the above example, but with $\frac{\partial u}{\partial y} = 0$ at the top boundary and $u = 0$ at the bottom boundary:
- -clear; clc; close all
-
-h = 0.1;
-x = [0 : h : 1]; n = length(x);
-y = [0 : h : 1]; m = length(y);
-
-% The coefficient matrix A is now m*n by m*n, since that is the total number of points.
-% The right-hand side vector b is m*n by 1.
-A = zeros(m*n, m*n);
-b = zeros(m*n, 1);
-
-u_left = 100;
-u_right = 100;
-u_bottom = 0;
-
-for j = 1 : m
- for i = 1 : n
- % for convenience we calculate all the indices once
- kij = (j-1)*n + i;
- kim1j = (j-1)*n + i - 1;
- kip1j = (j-1)*n + i + 1;
- kijm1 = (j-2)*n + i;
- kijp1 = j*n + i;
-
- if i == 1
- % this is the left boundary
- A(kij, kij) = 1;
- b(kij) = u_left;
- elseif i == n
- % right boundary
- A(kij, kij) = 1;
- b(kij) = u_right;
- elseif j == 1
- % bottom boundary
- A(kij, kij) = 1;
- b(kij) = u_bottom;
- elseif j == m
- % top boundary, using the ghost node + recursion formula
- A(kij, kim1j) = 1;
- A(kij, kip1j) = 1;
- A(kij, kijm1) = 2;
- A(kij, kij) = -4;
- else
- % these are the coefficients for the interior points,
- % based on the recursion formula
- A(kij, kim1j) = 1;
- A(kij, kip1j) = 1;
- A(kij, kijm1) = 1;
- A(kij, kijp1) = 1;
- A(kij, kij) = -4;
- end
- end
-end
-u = A \ b;
-
-u_square = reshape(u, [n, m]);
-% the "20" indicates the number of levels for the contour plot
-contourf(x, y, u_square', 20);
-c = colorbar;
-c.Label.String = 'Temperature';
-
-So far, we've been able to solve our systems of linear equations in Matlab by using y = A \ b, which directly finds the solution to the equation $A \mathbf{y} = \mathbf{b}$.
However, this approach will become very slow as the grid resolution ($h = \Delta x = \Delta y$) becomes smaller, and eventually unfeasable due to the associated computational requirements. First, let's create a function that takes as input the segment size $h$, then returns the time ittakes to solve the problem for different sizes.
- -%%file heat_equation.m
-function [time, num] = heat_equation(h)
-
-x = [0 : h : 1]; n = length(x);
-y = [0 : h : 1]; m = length(y);
-
-% The coefficient matrix A is now m*n by m*n, since that is the total number of points.
-% The right-hand side vector b is m*n by 1.
-A = zeros(m*n, m*n);
-b = zeros(m*n, 1);
-
-num = m*n; % number of points
-
-tic;
-
-u_left = 100;
-u_right = 100;
-u_bottom = 100;
-u_top = 0;
-
-for j = 1 : m
- for i = 1 : n
- % for convenience we calculate all the indices once
- kij = (j-1)*n + i;
- kim1j = (j-1)*n + i - 1;
- kip1j = (j-1)*n + i + 1;
- kijm1 = (j-2)*n + i;
- kijp1 = j*n + i;
-
- if i == 1
- % this is the left boundary
- A(kij, kij) = 1;
- b(kij) = u_left;
- elseif i == n
- % right boundary
- A(kij, kij) = 1;
- b(kij) = u_right;
- elseif j == 1
- % bottom boundary
- A(kij, kij) = 1;
- b(kij) = u_bottom;
- elseif j == m
- % top boundary
- A(kij, kij) = 1;
- b(kij) = u_top;
- else
- % these are the coefficients for the interior points,
- % based on the recursion formula
- A(kij, kim1j) = 1;
- A(kij, kip1j) = 1;
- A(kij, kijm1) = 1;
- A(kij, kijp1) = 1;
- A(kij, kij) = -4;
- end
- end
-end
-u = A \ b;
-
-u_square = reshape(u, [n, m]);
-% the "20" indicates the number of levels for the contour plot
-%contourf(x, y, u_square', 20);
-%c = colorbar;
-%c.Label.String = 'Temperature';
-
-time = toc;
-Now, we can see how long it takes to solve as we increase the resolution, and get an idea about the relationship between time-to-solution and number of unknowns.
- -clear all; clc;
-
-step_sizes = [0.1, 0.05, 0.025, 0.02, 0.0125, 0.01];
-
-n = length(step_sizes);
-nums = zeros(n,1); times = zeros(n,1);
-
-for i = 1 : n
- [times(i), nums(i)] = heat_equation(step_sizes(i));
-end
-
-loglog(nums, times, '-o')
-xlabel('Number of unknowns');
-ylabel('Time for direct solution (sec)')
-hold on
-
-x = nums(3:end);
-n2 = x.^2 * (times(3) / x(1)^2);
-n3 = x.^3 * (times(3) / x(1)^3);
-plot(x, n2, '-x');
-plot(x, n3, '-s');
-legend('Actual cost', 'Quadratic cost', 'Cubic cost', 'Location', 'northwest')
-
-Interestingly, we see that the slope in this log-log plot that after about 400 unknowns (so a coefficient matrix of about 160,000), the cost begins to increase exponentially, somewhere between quadratic ($\mathcal{O}(n^2)$) and cubic ($\mathcal{O}(n^3)$).
-If we try to reduce the step size further, for example to 0.005, we'll see that we cannot get a solution in a reasonable amount of time. But, clearly we want to get solutions for large numbers of unknowns, so what can we do?
-We can solve larger systems of linear equations using iterative methods. There are a number of these, and we'll focus on two:
-The Jacobi method essentially works by starting with an initial guess to the solution, then using the recursion formula to solve for values at each point, then repeating this until the values converge (i.e., stop changing).
-An algorithm we can use to solve Laplace's equation:
-More formally, if we have a system $A \mathbf{x} = \mathbf{b}$, where -\begin{equation} -A = \begin{bmatrix} -a_{11} & a_{12} & \cdots & a_{1n} \\ -a_{21} & a_{22} & \cdots & a_{2n} \\ -\vdots & \vdots & \ddots & \vdots \\ -a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} -\quad \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} -\quad \mathbf{b} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix} -\end{equation} -then we can solve iterative for $\mathbf{x}$ using -\begin{equation} -x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \right) , \quad i = 1,2,\ldots, n -\end{equation} -where $x_i^{(k)}$ is a value of the solution at iteration $k$ and $x_i^{(k+1)}$ is at the next iteration.
- -%%file heat_equation_jacobi.m
-function [time, num_point, num_iter] = heat_equation_jacobi(h)
-
-x = [0 : h : 1]; n = length(x);
-y = [0 : h : 1]; m = length(y);
-
-% The coefficient matrix A is now m*n by m*n, since that is the total number of points.
-% The right-hand side vector b is m*n by 1.
-A = zeros(m*n, m*n);
-b = zeros(m*n, 1);
-num_point = m*n;
-
-tic;
-
-u_left = 100;
-u_right = 100;
-u_bottom = 100;
-u_top = 0;
-
-% initial guess
-u = 100*ones(m*n, 1);
-
-% dummy value for residual variable
-epsilon = 1.0;
-
-num_iter = 0;
-while epsilon > 1e-6
- u_old = u;
-
- epsilon = 0;
- for j = 1 : m
- for i = 1 : n
- kij = (j-1)*n + i;
- kim1j = (j-1)*n + i - 1;
- kip1j = (j-1)*n + i + 1;
- kijm1 = (j-2)*n + i;
- kijp1 = j*n + i;
-
- if i == 1
- % this is the left boundary
- u(kij) = u_left;
- elseif i == n
- % right boundary
- u(kij) = u_right;
- elseif j == 1
- % bottom boundary
- u(kij) = u_bottom;
- elseif j == m
- % top boundary
- u(kij) = u_top;
- else
- % interior points
- u(kij) = (u_old(kip1j) + u_old(kim1j) + u_old(kijm1) + u_old(kijp1))/4.0;
- end
- end
- end
-
- epsilon = max(abs(u - u_old));
- num_iter = num_iter + 1;
-end
-
-u_square = reshape(u, [n, m]);
-% the "20" indicates the number of levels for the contour plot
-contourf(x, y, u_square', 20);
-c = colorbar;
-c.Label.String = 'Temperature';
-
-time = toc;
-fprintf('Number of iterations: %d\n', num_iter)
-step_sizes = [0.1, 0.05, 0.025, 0.02, 0.0125, 0.01, 0.005];
-n = length(step_sizes);
-
-nums_jac = zeros(n,1); times_jac = zeros(n,1); num_iter_jac = zeros(n,1);
-
-for i = 1 : n
- [times_jac(i), nums_jac(i), num_iter_jac(i)] = heat_equation_jacobi(step_sizes(i));
-end
-
-The Gauss-Seidel method is very similar to the Jacobi method, but with one important difference: rather than using all old values to calculate the new values, incorporate updated values as they are available. Because the method incorporates newer information more quickly, it tends to converge faster (meaning, with fewer iterations) than the Jacobi method.
-Formally, if we have a system $A \mathbf{x} = \mathbf{b}$, where -\begin{equation} -A = \begin{bmatrix} -a_{11} & a_{12} & \cdots & a_{1n} \\ -a_{21} & a_{22} & \cdots & a_{2n} \\ -\vdots & \vdots & \ddots & \vdots \\ -a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} -\quad \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} -\quad \mathbf{b} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix} -\end{equation} -then we can solve iterative for $\mathbf{x}$ using -\begin{equation} -x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j =i+1}^n a_{ij} x_j^{(k)} \right) , \quad i = 1,2,\ldots, n -\end{equation} -where $x_i^{(k)}$ is a value of the solution at iteration $k$ and $x_i^{(k+1)}$ is at the next iteration.
- -%%file heat_equation_gaussseidel.m
-function [time, num_point, num_iter] = heat_equation_gaussseidel(h)
-
-x = [0 : h : 1]; n = length(x);
-y = [0 : h : 1]; m = length(y);
-
-% The coefficient matrix A is now m*n by m*n, since that is the total number of points.
-% The right-hand side vector b is m*n by 1.
-A = zeros(m*n, m*n);
-b = zeros(m*n, 1);
-num_point = m*n;
-
-tic;
-
-u_left = 100;
-u_right = 100;
-u_bottom = 100;
-u_top = 0;
-
-% initial guess
-u = 100*ones(m*n, 1);
-
-% dummy value for residual variable
-epsilon = 1.0;
-
-num_iter = 0;
-while epsilon > 1e-6
- u_old = u;
-
- epsilon = 0;
- for j = 1 : m
- for i = 1 : n
- kij = (j-1)*n + i;
- kim1j = (j-1)*n + i - 1;
- kip1j = (j-1)*n + i + 1;
- kijm1 = (j-2)*n + i;
- kijp1 = j*n + i;
-
- if i == 1
- % this is the left boundary
- u(kij) = u_left;
- elseif i == n
- % right boundary
- u(kij) = u_right;
- elseif j == 1
- % bottom boundary
- u(kij) = u_bottom;
- elseif j == m
- % top boundary
- u(kij) = u_top;
- else
- % interior points
- u(kij) = (u(kip1j) + u(kim1j) + u(kijm1) + u(kijp1))/4.0;
- end
- end
- end
-
- epsilon = max(abs(u - u_old));
- num_iter = num_iter + 1;
-end
-
-u_square = reshape(u, [n, m]);
-%% the "20" indicates the number of levels for the contour plot
-contourf(x, y, u_square', 20);
-c = colorbar;
-c.Label.String = 'Temperature';
-
-time = toc;
-fprintf('Number of iterations: %d\n', num_iter)
-step_sizes = [0.1, 0.05, 0.025, 0.02, 0.0125, 0.01, 0.005];
-n = length(step_sizes);
-
-nums_gs = zeros(n,1); times_gs = zeros(n,1); num_iter_gs = zeros(n,1);
-
-for i = 1 : n
- [times_gs(i), nums_gs(i), num_iter_gs(i)] = heat_equation_gaussseidel(step_sizes(i));
-end
-
-loglog(nums, times, '-o'); hold on
-loglog(nums_jac, times_jac, '-^')
-loglog(nums_gs, times_gs, '-x')
-xlabel('Number of unknowns');
-ylabel('Time for solution (sec)')
-legend('Direct solution', 'Jacobi solution', 'Gauss-Seidel solution', 'Location', 'northwest')
-
-loglog(nums_jac, num_iter_jac, '-^')
-hold on;
-loglog(nums_gs, num_iter_gs, '-x')
-xlabel('Number of unknowns')
-ylabel('Number of iterations required')
-legend('Jacobi method', 'Gauss-Seidel method', 'Location', 'northwest')
-
-These results show us a few things:
-For larger, more-realistic problems, iterative solution methods like Jacobi and Gauss-Seidel are essential.
- -A classic example of a parabolic partial differential equation (PDE) is the one-dimensional unsteady heat equation: -\begin{equation} -\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial t^2} -\end{equation} -where $T(x,t)$ is the temperature varying in space and time, and $\alpha$ is the thermal diffusivity: $\alpha = k / (\rho c_p)$, which is a constant.
-We can solve this using finite differences to represent the spatial derivatives and time derivatives separately. -First, let's rearrange the PDE slightly: -\begin{equation} -\frac{\partial^2 T}{\partial x^2} = \frac{1}{\alpha} \frac{\partial T}{\partial t} -\end{equation}
-Let's use a central difference for the spatial derivative with a spacing of $\Delta x$, and a forward difference for the time derivative with a time-step size of $\Delta t$. With these choices, we can obtain an approximation to the PDE that applies at time $t^k$ and spatial location $x_i$: -\begin{equation} -\frac{T_{i-1}^k - 2 T_i^k + T_{i+1}^k}{\Delta x^2} = \frac{1}{\alpha} \left( \frac{T_i^{k+1} - T_i^k}{\Delta t} \right) -\end{equation} -where $T_i^k$ is the temperature at time $t^k$ and spatial location $x_i$. The following figure shows the stencil of points involved in the PDE, for a domain with five points in the $x$-direction.
-
- To solve the heat equation for a one-dimensional domain over $0 \leq x \leq L$, we will need both initial conditions at $t = 0$ and boundary conditions at $x=0$ and $x=L$ (for all time). In terms of our nodal values, this means we need $T_i^{k=1}$ for $i = 1 \ldots n$, where $n$ is the number of points, as well as information about $T_1^k$ and $T_n^k$ for all times $k$.
-Based on the choices of finite differences, we have some information on the accuracy of this approach:
-We can rearrange the above equation to obtain our recursion formula: -\begin{equation} -T_i^{k+1} = \left( T_{i+1}^k + T_{i-1}^k \right) \frac{\alpha \Delta t}{\Delta x^2} + T_i^k \left( 1 - 2 \frac{\alpha \Delta t}{\Delta x^2} \right) \;. -\end{equation} -This is an explicit scheme in time, similar to the Forward Euler method we used for ordinary differential equations, and like that method it may have stability issues. The combination of terms we see repeated is also known as the Fourier number: $\text{Fo} = \frac{\alpha \Delta t}{\Delta x^2}$, and governs the stability. -We can rewrite the recursion formula using this: -\begin{equation} -T_i^{k+1} = \left( T_{i+1}^k + T_{i-1}^k \right) \text{Fo} + T_i^k \left( 1 - 2 \text{Fo} \right) \;. -\end{equation}
-The term in parentheses there must be greater than or equal to zero for stability ($1 - 2 \text{Fo} \geq 0$); if not, the solution may become unstable and blow up. This gives us some conditions on our choice of step sizes: -\begin{align} -1 - 2 \text{Fo} &\geq 0 \\ -1 & \geq 2 \text{Fo} \\ -\text{Fo} &\leq \frac{1}{2} \\ -\frac{\alpha \Delta t}{\Delta x^2} &\leq \frac{1}{2} -\end{align} -This is the stability criterion for the explicit method: the Fourier number must be smaller than 0.5. -For a given thermal diffusivity and chosen spatial step size, this also gives us the limit on the time-step size: $\Delta t \leq \Delta x^2 / (2 \alpha)$.
-Let's look at an example where the initial temperature is 200, the temperature at the boundaries are 50, the thermal diffusivity is $\alpha = 2.3 \times 10^{-1}$ m$^2 /$ s, and $L = 1$.
-In other words,
-\begin{align}
-T(x, t=0) &= 200 \\
-T(x=0, t) &= 50 \\
-T(x=L, t) &= 50
-\end{align}
-We'll integrate out to $t = 1$, using a Fourier number of 0.25 to be comfortably below the stability limit (Fo = 0.25):
clear all
-
-dx = 0.1;
-alpha = 2.3e-1;
-
-% for stability, set the Fourier number at 0.25 (half the stability limit of 0.5)
-Fo = 0.25;
-% then choose the time step based on the Fourier number
-dt = Fo * dx^2 / alpha;
-
-fprintf('Time step size: %5.3f\n', dt);
-
-x = [0 : dx : 1]; n = length(x);
-t = [0 : dt : 1]; m = length(t);
-
-T = zeros(m, n);
-
-% initial conditions
-T(1,:) = 200;
-
-plot(x, T(1,:))
-axis([0 1 50 200]);
-xlabel('Distance'); ylabel('Temperature');
-F(1) = getframe(gcf);
-
-for k = 1 : m - 1
- for i = 1 : n
- if i == 1
- T(k+1, 1) = 50;
- elseif i == n
- T(k+1, n) = 50;
- else
- T(k+1, i) = (T(k,i+1) + T(k,i-1))*Fo + T(k,i)*(1 - 2*Fo);
- end
- end
- plot(x, T(k+1,:))
- axis([0 1 50 200]);
- xlabel('Distance'); ylabel('Temperature');
- F(k+1) = getframe(gcf);
-end
-close
-
-%% If you are working interactively, you can use this to make a movie in Matlab
-%fig = figure;
-%movie(fig, F, 2)
-
-%% This generates a GIF of the results (for use in Jupyter Notebook)
-filename = 'parabolic_animated.gif';
-for i = 1 : length(F)
- im = frame2im(F(i));
- [imind,cm] = rgb2ind(im,256);
- % Write to GIF
- if i == 1
- imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3);
- else
- imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
- end
-end
-
- This shows the temperature decaying exponentially from the initial conditions, constrained by the boundary conditions.
-What happens if we tried to use a Fourier number larger than 0.5, or arbitrarily chose a time-step size that was too large (and resulted in $\text{Fo} > 0.5$)?
- -clear all
-
-dx = 0.1;
-alpha = 2.3e-1;
-
-%% Purposely choose a Fourier number that is past the stability limit:
-Fo = 0.75;
-dt = Fo * dx^2 / alpha;
-
-fprintf('Time step size: %5.3f\n', dt);
-
-x = [0 : dx : 1]; n = length(x);
-t = [0 : dt : 1]; m = length(t);
-
-T = zeros(m, n);
-
-% initial conditions
-T(1,:) = 200;
-
-plot(x, T(1,:))
-axis([0 1 50 200]);
-xlabel('Distance'); ylabel('Temperature');
-F(1) = getframe(gcf);
-
-for k = 1 : m - 1
- for i = 1 : n
- if i == 1
- T(k+1, 1) = 50;
- elseif i == n
- T(k+1, n) = 50;
- else
- T(k+1, i) = (T(k,i+1) + T(k,i-1))*Fo + T(k,i)*(1 - 2*Fo);
- end
- end
- plot(x, T(k+1,:))
- xlabel('Distance'); ylabel('Temperature');
- F(k+1) = getframe(gcf);
-end
-close
-
-%% If you are working interactively, you can use this to make a movie in Matlab
-%fig = figure;
-%movie(fig, F, 2)
-
-%% This generates a GIF of the results (for use in Jupyter Notebook)
-filename = 'parabolic_unstable_animated.gif';
-for i = 1 : length(F)
- im = frame2im(F(i));
- [imind,cm] = rgb2ind(im,256);
- % Write to GIF
- if i == 1
- imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3);
- else
- imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
- end
-end
-
- In this case, the solution becomes unstable and blows up, leading to unphysical results.
-For this explicit scheme, which is conditionally stable, the choice of $\Delta t$ is limited by the stability criterion. This means that we may be stuck using a small time-step size.
-Rather than being forced to use a very small time-step size, we can also explore implicit schemes that are unconditionally stable.
- -To solve the problem without being constrained by the stability criterion (and the associated restriction on time-step size), we can develop an implicit scheme for solving the 1D unsteady heat equation by evaluating the second derivative at time $t^{k+1}$, rather than at time $t^k$ as we did above: -\begin{equation} -\frac{T_i^{k+1} - T_i^k}{\Delta t} = \alpha \frac{T_{i-1}^{k+1} - 2 T_i^{k+1} + T_{i+1}^{k+1}}{\Delta x^2} -\end{equation}
-Then, rearrange this to obtain a recursion formula with unknowns on the left-hand side and knowns on the right-hand side: -\begin{equation} -T_{i-1}^{k+1} (\text{Fo}) + T_i^{k+1} \left(-2 \text{Fo} - 1 \right) + T_{i+1}^{k+1} (\text{Fo}) = -T_i^k -\end{equation} -where $\text{Fo} = \alpha \Delta t / \Delta x^2$ is the Fourier number.
-
- Unlike the explicit scheme, with the implicit scheme we do not have a simple recursion formula that can be applied to calculate the next temperature value at each location. Instead, we need to solve a system of linear equations at each time step, to simultaneously get all the values of temperature at the next time step. -In other words, at each time step, we need to solve -\begin{equation} -A \mathbf{T}^{k+1} = - \mathbf{T}^k = \mathbf{b} \;. -\end{equation}
-Let's implement this to solve the same example:
- -clear all
-
-alpha = 2.3e-1;
-dx = 0.1;
-x = [0 : dx : 1]; n = length(x);
-
-% choose a Fourier number that is deliberately past the explicit method stability limit
-Fo = 0.75;
-dt = Fo * dx^2 / alpha;
-
-fprintf('Time step size: %5.3f\n', dt);
-
-t = [0 : dt : 1]; m = length(t);
-
-T = zeros(m, n);
-
-% Initial conditions
-T(1,:) = 200;
-
-plot(x, T(1,:))
-axis([0 1 50 200]);
-xlabel('Distance'); ylabel('Temperature');
-F(1) = getframe(gcf);
-
-for k = 1 : m - 1
- A = zeros(n,n);
- b = zeros(n,1);
- for i = 1 : n
- if i == 1
- A(1,1) = 1;
- b(1) = 50;
- elseif i == n
- A(n,n) = 1;
- b(n) = 50;
- else
- A(i,i-1) = Fo;
- A(i,i) = -2*Fo - 1;
- A(i,i+1) = Fo;
- b(i) = -T(k,i);
- end
- end
-
- T(k+1, :) = A \ b;
- plot(x, T(k+1,:))
- axis([0 1 50 200]);
- xlabel('Distance'); ylabel('Temperature');
- F(k+1) = getframe(gcf);
-end
-close
-
-%% If you are working interactively, you can use this to make a movie in Matlab
-%fig = figure;
-%movie(fig, F, 2)
-
-%% This generates a GIF of the results (for use in Jupyter Notebook)
-filename = 'parabolic_implicit_animated.gif';
-for i = 1 : length(F)
- im = frame2im(F(i));
- [imind,cm] = rgb2ind(im,256);
- % Write to GIF
- if i == 1
- imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3);
- else
- imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
- end
-end
-
- With this implicit scheme, we can now use larger time-step sizes (resulting in values of the Fourier number greater than 0.5), because the method is unconditionally stable. Like the explicit scheme above, this approach is:
-So far we have two methods to solve parabolic PDEs:
-We might also want a scheme that is more accurate with respect to the time-step size $\Delta t$. The Crank-Nicolson method is 2nd-order accurate in time and also unconditionally stable (because it is implicit). This method is essentially an average of the first-order accurate explicit and implicit methods: -\begin{equation} -\frac{T_i^{k+1} - T_i^k}{\Delta t} = \frac{\alpha}{2} \frac{T_{i-1}^{k} - 2 T_i^{k} + T_{i+1}^{k}}{\Delta x^2} + \frac{\alpha}{2} \frac{T_{i-1}^{k+1} - 2 T_i^{k+1} + T_{i+1}^{k+1}}{\Delta x^2} -\end{equation} -This method is based on using the average of the time derivative at the current and next time steps.
-
- If we rearrange that equation, we can get our recursion formula: -\begin{equation} -T_{i-1}^{k+1} \left( \frac{\text{Fo}}{2} \right) + T_{i}^{k+1} \left( Fo - 1 \right) + T_{i+1}^{k+1} \left( \frac{\text{Fo}}{2} \right) = T_{i-1}^{k} \left( -\frac{\text{Fo}}{2} \right) + T_{i}^{k} \left( \text{Fo} - 1 \right) + T_{i-1}^{k} \left( -\frac{\text{Fo}}{2} \right) \;, -\end{equation} -where $\text{Fo} = \alpha \Delta t / \Delta x^2$. -Like the first-order implicit scheme, advancing this solution in time requires solving a system of linear equations at each time step: -\begin{equation} -A \mathbf{T}^{k+1} = B \mathbf{T}^k = \mathbf{b} \;, -\end{equation} -where the coefficient matrices $A$ and $B$ will look like -\begin{equation} -A = \begin{bmatrix} -1 \\ -\frac{\text{Fo}}{2} & (-\text{Fo}-1) & \frac{\text{Fo}}{2} \\ - & \frac{\text{Fo}}{2} & (-\text{Fo}-1) & \frac{\text{Fo}}{2} \\ - & & \ddots & \ddots & \ddots & \\ - & & & & & 1 -\end{bmatrix} -\quad -B = \begin{bmatrix} -1 \\ -\frac{-\text{Fo}}{2} & (\text{Fo}-1) & \frac{\text{-Fo}}{2} \\ - & \frac{-\text{Fo}}{2} & (\text{Fo}-1) & \frac{-\text{Fo}}{2} \\ - & & \ddots & \ddots & \ddots & \\ - & & & & & 1 -\end{bmatrix} -\end{equation} -The first and last rows of these matrices may differ, depending on the boundary conditions.
-Let's implement this method:
- -clear all
-
-alpha = 2.3e-1;
-dx = 0.1;
-x = [0 : dx : 1]; n = length(x);
-
-% choose a Fourier number that is deliberately past the explicit method stability limit
-Fo = 0.75;
-dt = Fo * dx^2 / alpha;
-
-fprintf('Time step size: %5.3f\n', dt);
-
-t = [0 : dt : 1]; m = length(t);
-
-T = zeros(m, n);
-
-% Initial conditions
-T(1,:) = 200;
-
-% boundary conditions
-T(:,1) = 50;
-T(:,n) = 50;
-
-plot(x, T(1,:))
-axis([0 1 50 200]);
-xlabel('Distance'); ylabel('Temperature');
-F(1) = getframe(gcf);
-
-for k = 1 : m - 1
- A = zeros(n,n);
- B = zeros(n,n);
- for i = 1 : n
- if i == 1
- A(1,1) = 1;
- B(1,1) = 1;
- elseif i == n
- A(n,n) = 1;
- B(n,n) = 1;
- else
- A(i,i-1) = Fo/2;
- A(i,i) = -Fo - 1;
- A(i,i+1) = Fo/2;
- B(i,i-1) = -Fo/2;
- B(i,i) = Fo - 1;
- B(i,i+1) = -Fo/2;
- end
- end
- b = B * T(k,:)';
- T(k+1, :) = A \ b;
- plot(x, T(k+1,:))
- axis([0 1 50 200]);
- xlabel('Distance'); ylabel('Temperature');
- F(k+1) = getframe(gcf);
-end
-close
-
-%% If you are working interactively, you can use this to make a movie in Matlab
-%fig = figure;
-%movie(fig, F, 2)
-
-%% This generates a GIF of the results (for use in Jupyter Notebook)
-filename = 'parabolic_cranknicolson_animated.gif';
-for i = 1 : length(F)
- im = frame2im(F(i));
- [imind,cm] = rgb2ind(im,256);
- % Write to GIF
- if i == 1
- imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3);
- else
- imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
- end
-end
-
- This chapter focuses on numerical methods for solving partial differential equations (PDEs), which involve derivatives in multiple dimensions.
-We can write a general, linear 2nd-order PDE for a variable $u(x,y)$ as -\begin{equation} -A \frac{\partial^2 u}{\partial x^2} + 2 B \frac{\partial^2 u}{\partial x \, \partial y} + C \frac{\partial^2 u}{\partial y^2} = F \left( x, y, u, \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y} \right) -\end{equation} -where $A$, $B$, and $C$ are constants. Depending on their value, we can categorize a PDE into one of three categories:
-The different PDE types will exhibit different characteristics and will also require slightly different solution approaches.
- -Given the first-order ODE -\begin{equation} -\frac{dy}{dt} = -8y \;, -\end{equation}
-a.) Perform a linear stability analysis to find the amplification factor ($\sigma$) for the midpoint method (also known as the modified Euler method) applied to this ODE.
-b.) Using your result from part (a), show whether the numerical solution would be stable or unstable for time-step sizes of $\Delta t = 0.2$ and $\Delta t = 0.4$. (Show for each value)
-c.) Based on your results for parts (a) and (b), is the midpoint method unstable, conditionally stable, or unconditionally stable for this ODE?
-d.) What is the order of accuracy for the midpoint method? Based on this, what (approximate) global errors would you expect in the solution when using time-step sizes of $\Delta t = 0.2$ and 0.4?
-e.) What are your two options for reducing error in the solution?
-a.) -\begin{align} -y_{i+1/2} &= y_i + \frac{\Delta t}{2} \, f(t_i, y_i) \\ -y_{i+1} &= y_i + \Delta t \, f(t_i + \frac{\Delta t}{2}, y_{i+1/2}) -\end{align}
-For this ODE: -\begin{align} -y_{i+1/2} &= y_i + \frac{\Delta t}{2} \, (-8 y_i) = y_i - 4 \Delta t \, y_i \\ -y_{i+1} &= y_i + \Delta t \, \left[ -8 \left( y_i - 4 \Delta t \, y_i \right) \right] \\ -&= y_i - 8 \Delta t \, y_i + 32 \Delta t^2 \, y_i \\ -&= y_i (1 - 8 \Delta t + 32 \Delta t^2) -\end{align} -So, -\begin{equation} -\sigma = \frac{y_{i+1}}{y_i} = 1 - 8 \Delta t + 32 \Delta t^2 -\end{equation}
-b.) For stability, $|\sigma| \leq 1$.
-$\Delta t = 0.2$: -\begin{equation} -\sigma = 1 - 8(0.2) + 32 (0.2)^2 = 0.68 -\end{equation} -so stable.
-$\Delta t = 0.4$: -\begin{equation} -\sigma = 1 - 8(0.4) + 32 (0.4)^2 = 2.92 -\end{equation} -so unstable.
-c.) Conditionally stable; the method is stable for some values of $\Delta t$ and unstable for other values.
-d.) The midpoint method is 2nd-order accurate. So, for $\Delta t = 0.2$, we should expect global errors on the order of 0.04, and for $\Delta t = 0.4$ we should expect errors on the order of 0.16.
-e.)
-Given the second-order ODE -\begin{equation} -2 y^{\prime\prime} + y^{\prime} + 4y = 3x -\end{equation}
-a.) Find the recursion formulas (i.e., $y_{i+1} = \ldots$) for numerically solving this ODE using the backward Euler method. Clearly define/state any variable or function you use.
-b.) What is the order of accuracy for the backward Euler method? Given a step size $\Delta x = 0.15$, approximately what local error and what global error would you expect in your solution? What is the difference between these two errors?
-c.) Why would you want to use this method to solve an ODE over a simpler method like forward Euler?
-a.) System of Backward Euler recursion formulas, for $y(x)$ and $u(x) = y^{\prime}$, where $f(x,y,u) = dy/dx$ and $g(x,y,u) = du/dx = y^{\prime\prime}$: -\begin{align} -y_{i+1} &= y_i + \Delta x \, f(x_{i+1}, y_{i+1}, u_{i+1}) = y_i + \Delta x \, u_{i+1} \\ -u_{i+1} &= u_i + \Delta x \, g(x_{i+1}, y_{i+1}, u_{i+1}) = u_i + \Delta x \, \left( \frac{3}{2}x_{i+1} - \frac{1}{2} u_{i+1} - 2 y_{i+1} \right) -\end{align} -This form is implicit and cannot be used directly, so we need to solve the system of equations.
-There are two ways to approach this, that give an equivalent solution
-Option 1: Use substitution and elimination to solve: -\begin{align} -u_{i+1} &= u_i + \frac{3}{2} \Delta x \, x_{i+1} - \frac{1}{2} \Delta x \, u_{i+1} - 2 \Delta x \left( y_i + \Delta x \, u_{i+1} \right) \\ -&= u_i + \frac{3}{2} \Delta x \, x_{i+1} - \frac{1}{2} \Delta x \, u_{i+1} - 2 \Delta x \, y_i - 2 \Delta x^2 \, u_{i+1} \\ -u_{i+1} \left( 1 + 2 \Delta x^2 + \frac{1}{2} \Delta x \right) &= u_i + \frac{3}{2} \Delta x \, x_{i+1} - 2 \Delta x \, y_i -\end{align} -Thus, -\begin{align} -u_{i+1} &= \frac{u_i + \frac{3}{2} \Delta x \, x_{i+1} - 2 \Delta x \, y_i}{1 + 2 \Delta x^2 + \frac{1}{2} \Delta x} \\ -y_{i+1} &= y_i + \Delta x \frac{u_i + \frac{3}{2} \Delta x \, x_{i+1} - 2 \Delta x \, y_i}{1 + 2 \Delta x^2 + \frac{1}{2} \Delta x} -\end{align}
-Option 2: Or, use Cramer's rule: -\begin{align} -y_{i+1} - \Delta x u_{i+1} &= y_i \\ -2 \Delta x \, y_{i+1} + \left( 1 + \frac{1}{2} \Delta x \right) u_{i+1} &= u_i + \frac{3}{2} \Delta x \, x_{i+1} \\ -\rightarrow \begin{bmatrix} 1 & -\Delta x \\ 2 \Delta x & \left( 1 + \frac{1}{2} \Delta x \right)\end{bmatrix} \begin{bmatrix} y_{i+1} \\ u_{i+1} \end{bmatrix} &= -\begin{bmatrix} y_i \\ u_i + \frac{3}{2} \Delta x \, x_{i+1} \end{bmatrix} -\end{align} -Then, -\begin{align} -y_{i+1} &= \frac{y_i \left(1 + \frac{\Delta x}{2} \right) + \Delta x \left( u_i + \frac{3}{2} \Delta x \, x_{i+1} \right)}{1 + \frac{\Delta x}{2} + 2 \Delta x^2} \\ -u_{i+1} &= \frac{u_i + \frac{3}{2} \Delta x \, x_{i+1} - 2 \Delta x \, y_i}{1 + \frac{\Delta x}{2} + 2 \Delta x^2} -\end{align}
-b.) Backward Euler is 1st-order accurate, so the global error is on the order of the step size ($\Delta x$).
-The global error should then be on the order of 0.15, and the local error on the order of $0.15^2 = 0.0225$.
-The difference: the local (or truncation) error is the error at each step of the solution, while the global error is the overall error that accumulates over the whole solution.
-c.) Backward Euler is unconditionally stable, while Forward Euler is conditionally stable.
-Given the input waveform $R(t)$ shown here,
-
a.) What is the period and fundamental frequency of the input forcing function?
-b.) Is the periodic function $R(t)$ odd, even, or neither?
-c.) Find the coefficients $a_0$ and $a_n$ of the Fourier series representation of $R(t)$. (For purposes of time, you do not need to find $b_n$): -\begin{equation} -R(t) = a_0 + \sum_{n=1}^{\infty} a_n \cos(n \omega t) + \sum_{n=1}^{\infty} b_n \sin(n \omega t) -\end{equation}
-a.) -\begin{equation} -T = 4\pi \quad \omega = \frac{2\pi}{T} = \frac{1}{2} -\end{equation}
-b.) Neither.
-c.) -\begin{align} -a_0 &= \frac{1}{T} \int_0^T f(t) dt = \frac{1}{4\pi} \left[ \int_0^{\pi} 0 dt + \int_{\pi}^{2\pi} 3 dt + \int_{2\pi}^{3\pi} 2 dt + \int_{3\pi}^{4\pi} 1 dt \right] \\ -&= \frac{1}{4\pi} \left[ 3 \pi + 2 \pi + 1 \pi \right] \\ -a_0 &= \frac{3}{2} -\end{align}
-\begin{align} -a_n &= \frac{2}{T} \int_0^T f(t) \cos(n \omega t) dt \\ -&= \frac{2}{4\pi} \left[ \int_0^{\pi} 0 \cos \left(\frac{n t}{2}\right) dt + \int_{\pi}^{2\pi} 3 \cos \left(\frac{n t}{2}\right) dt + \int_{2\pi}^{3\pi} 2 \cos \left(\frac{n t}{2}\right) dt + \int_{3\pi}^{4\pi} 1 \cos \left(\frac{n t}{2}\right) dt \right] \\ -&\cdots \\ -a_n &= \frac{1}{n\pi} \left[ \sin\left(\frac{3 n \pi}{2}\right) - 3 \sin\left(\frac{n \pi}{2}\right) \right] -\end{align}The displacement $y(t)$ of a harmonically forced mass-spring system is given by: -\begin{equation} -y^{\prime\prime} + 8y^{\prime} + 16y = 6 e^{-4t} -\end{equation}
-a.) For initial conditions $y(0)=0$ and $y^{\prime}(0) = 2$, find the response of the system $y(t)$.
-b.) Given a specified time increment $\Delta t$ and a domain, write the recursion formulas for solving this equation with the forward Euler method. Clearly define any variables or functions used. You do not need to write Matlab code.
-a.) First, get the homogeneous solution: -\begin{align} -y_H^{\prime\prime} + 8y_H^{\prime} + 16y_H &= 0 \\ -\rightarrow \lambda^2 + 8 \lambda + 16 &= 0 = (\lambda + 4)^2 \\ -\lambda &= 4 \text{ (repeated)} \\ -\text{so } y_H &= c_1 e^{-4t} + c_2 t e^{-4t} -\end{align}
-Next, use the method of undetermined coefficients to get the inhomogeneous solution: -\begin{align} -y_{IH} &= K t^2 e^{-4t} \\ -y_{IH}^{\prime} &= K e^{-4t} \left( 2t - 4t^2 \right) \\ -y_{IH}^{\prime\prime} &= K e^{-4t} \left( -8t + 2 - 8t + 16t^2 \right) \\ -\rightarrow K &= 3 -\end{align}
-Then, the general solution is -\begin{equation} -y(t) = c_1 e^{-4t} + c_2 t e^{-4t} + 3 t^2 e^{-4t} -\end{equation}
-Applying the initial conditions: -\begin{equation} -y(t) = 2 t e^{-4t} + 3 t^2 e^{-4t} -\end{equation}
-b.) If $z_1 = y$ and $z_2 = y^{\prime}$, then -\begin{align} -z_1^{\prime} &= z_2 \\ -z_2^{\prime} &= y^{\prime\prime} = 6e^{-4t} - 8z_2 - 16z_1 = f(t, z_1, z_2) -\end{align}
-Then, the recursion formulas are: -\begin{align} -z_{1, i+1} &= z_{1,i} + \Delta t \, z_{2,i} \\ -z_{2, i+1} &= z_{2, i} + \Delta t \left( 6 e^{-4t} - 8 z_{2,i} - 16 z_{1,i} \right) -\end{align} -where $z_{1,1} = 0$ and $z_{2,1} = 4$.
- -The temperature distribution $T(r)$ in an annular fin of inner radius $r_1$ and outer radius $r_2$ is described by the equation -\begin{equation} -r \frac{d^2 T}{dr^2} + \frac{dT}{dr} - rm^2 (T - T_{\infty}) = 0 \;, -\end{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: -\begin{equation} --k \left. \frac{dT}{dr} \right|_{r=r_2} = h \left[ T(r=r_2) - T_{\infty} \right] \;. -\end{equation} -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: -\begin{equation} -r_i \frac{T_{i-1} - 2T_i + T_{i+1}}{\Delta r^2} + \frac{T_{i+1} - T_{i-1}}{2\Delta r} - r_i m^2 (T_i - T_{\infty}) = 0 -\end{equation} -or -\begin{equation} -r_i (T_{i-1} - 2T_i + T_{i+1}) + \frac{\Delta r}{2} (T_{i+1} - T_{i-1}) - r_i m^2 \Delta r^2 (T_i - T_{\infty}) = 0 -\end{equation}
-b.) Rearrange and combine terms: -\begin{align} -r_i (T_{i-1} - 2T_i + T_{i+1}) + \frac{\Delta r}{2} (T_{i+1} - T_{i-1}) - r_i m^2 \Delta r^2 (T_i - T_{\infty}) &= 0 \\ -r_i (T_{i-1} - 2T_i + T_{i+1}) + \frac{\Delta r}{2} (T_{i+1} - T_{i-1}) - r_i m^2 \Delta r^2 T_i &= -r_i m^2 \Delta r^2 T_{\infty} \\ -\left(r_i - \frac{\Delta r}{2}\right) T_{i-1} + \left( -2 r_i - r_i m^2 \Delta r^2 \right) T_i + \left( r_i + \frac{\Delta r}{2} \right) T_{i+1} &= -r_i m^2 \Delta r^2 T_{\infty} -\end{align}
-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$: -\begin{align} --k \frac{T_n - T_{n-1}}{\Delta r} &= h (T_n - T_{\infty}) \\ --k (T_n - T_{n-1}) &= h \Delta r (T_n - T_{\infty}) \\ -k T_{n-1} - (k + h\Delta r) T_n &= -h \Delta r T_{\infty} -\end{align}
-Given the equation $y^{\prime\prime} + 9 \lambda^2 y = 0$ with $y(0) = 0$ and $y(2) = 0$,
-a.) Find the expression that gives all eigenvalues ($\lambda$). What is the eigenfunction?
-b.) Calculate the principal eigenvalue.
-a.) -\begin{gather} -y(x) = A \sin (3 \lambda x) + B \cos (3 \lambda x) \\ -\text{Apply BCs: } y(x=0) = 0 = A \sin(0) + B \cos(0) = B \\ -\therefore B = 0 \\ -y(x) = A \sin (3 \lambda x) \\ -y(x=2) = 0 = A \sin (3 \lambda 2) \\ -A \neq 0 \text{ so } \sin(3 \lambda 2) = \sin(6 \lambda) = 0 \therefore 6 \lambda = n \pi \quad n=1,2,3,\ldots \\ -\lambda = \frac{n \pi}{6} \quad n=1,2,3,\ldots,\infty -\end{gather}
-The eigenfunction is then the solution function associated with an eigenvalue: -\begin{equation} -y_n = A_n \sin \left( \frac{n \pi x}{2} \right) \quad n = 1, 2, 3, \ldots, \infty -\end{equation}
-b.) The principal eigenvalue is just that associated with $n = 1$: -\begin{equation} -\lambda_p = \lambda_1 = \frac{\pi}{6} -\end{equation}
-Use the shooting method to solve the boundary value problem -\begin{equation} -y^{\prime\prime} - 4y = 0 -\end{equation} -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: -\begin{align} -z_1' &= y' = z_2 \\ -z_2' &= y'' = 4 z_1 -\end{align} -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: -\begin{align} -z_{1,2} = z_1 (0.5) &= z_1 (0) + z_2(0) 0.5 = 0 \\ -z_{2,2} = z_2 (0.5) &= z_2 (0) + \left( 4z_1(0) \right) 0.5 = 0 \\ -z_{1,3} = z_1 (1.0) &= z_1 (0.5) + z_2(0.5) 0.5 = 0 \leftarrow \text{solution 1} \\ -z_{2,3} = z_2 (1.0) &= z_2 (0.5) + \left( 4z_1(0.5) \right) 0.5 = 0 -\end{align} -so for solution 1: $y(1) = 0 \neq 3$.
-For guess #2: $y' (0) = 2 = z_2 (0)$, with the forward Euler method: -\begin{align} -z_1 (0.5) &= z_1 (0) + z_2(0) 0.5 = 1.0 \\ -z_2 (0.5) &= z_2 (0) + \left( 4z_1(0) \right) 0.5 = 2.0 \\ -z_1 (1.0) &= z_1 (0.5) + z_2(0.5) 0.5 = 2.0 \leftarrow \text{solution 2} \\ -z_2 (1.0) &= z_2 (0.5) + \left( 4z_1(0.5) \right) 0.5 = 4.0 -\end{align} -so for solution 1: $y(1) = 2 \neq 3$.
-For guess #3, we can interpolate: -\begin{align} -m &= \frac{\text{guess 1} - \text{guess 2}}{\text{solution 1} - \text{solution 2}} = \frac{0 - 2}{0 - 2} = 1 \\ -\text{guess 3} &= \text{guess 2} + m (\text{target} - \text{solution 2}) = 2 + 1(3-2) = 3 -\end{align} -then, use this guess: -\begin{align} -z_1 (0.5) &= z_1 (0) + z_2(0) 0.5 = 1.5 \\ -z_2 (0.5) &= z_2 (0) + \left( 4z_1(0) \right) 0.5 = 3.0 \\ -z_1 (1.0) &= z_1 (0.5) + z_2(0.5) 0.5 = 3.0 \leftarrow \text{solution 3} \\ -z_2 (1.0) &= z_2 (0.5) + \left( 4z_1(0.5) \right) 0.5 = 6.0 -\end{align} -so for solution 3: $y(1) = 3$ which is the target.
-So our answer is $y'(0) = 3$.
-
-Let's first focus on simple analytical solutions for 2nd-order ODEs. As before, let's categorize problems based on their solution approach.
-If you have a 2nd-order ODE of this form: -\begin{equation} -\frac{d^2 y}{dx^2} = f(x) -\end{equation} -then you can solve by direct integration.
-For example, let's say we are trying to solve for the deflection of a cantilever beam $y(x)$ with a force $P$ at the end, where $E$ is the modulus and $I$ is the moment of intertia, and the initial conditions are $y(0)=0$ and $y^{\prime}(0) = 0$:
-
That is our general solution; we can obtain the specific solution by applying our two initial conditions: -\begin{align} -y(0) &= 0 = C_2 \\ -y^{\prime}(0) &= 0 = C_1 \\ -\therefore y(x) &= \frac{P}{EI} \left( \frac{x^3}{6} - \frac{L x^2}{2} \right) -\end{align}
- -If we have a 2nd-order ODE of this form: -\begin{equation} -\frac{d^2 y}{dx^2} = f(x, y^{\prime}) -\end{equation} -then we can solve by substitution, meaning by substituting a new variable for $y^{\prime}$. (Notice that $y$ itself does not show up in the ODE.)
-Let's substitute $u$ for $y^{\prime}$ in the above ODE: -\begin{align} -u &= y^{\prime} \\ -u^{\prime} &= y^{\prime\prime} \\ -\rightarrow u^{\prime} &= f(f, u) -\end{align} -Now we have a 1st-order ODE! Then, we can apply the methods previously discussed to solve this; once we find $u(x)$, we can integrate that once more to get $y(x)$.
-For example, consider a falling mass where we are solving for the downward distance as a function of time, $y(t)$, that is experiencing the force of gravity downward and a drag force upward. It starts at some reference point so $y(0) = 0$, and has a zero initial downward velocity: $y^{\prime}(0) = 0$. The governing equation is: -\begin{equation} -m \frac{d^2 y}{dt^2} = mg - c \left( \frac{dy}{dt} \right)^2 -\end{equation} -where $m$ is the mass, $g$ is acceleration due to gravity, and $c$ is the drag proportionality constant. -We can substitute $V$ for $y^{\prime}$, which gives us a first-order ODE: -\begin{align} -\text{let} \quad \frac{dy}{dt} = V \\ -\rightarrow m \frac{dV}{dt} &= mg - c V^2 \\ -\end{align} -Then, we can solve this for $V(t)$ using our initial condition for velocity $V(0) = 0$. Once we have that, we can integrate once more: -\begin{equation} -y(t) = \int V(t) dt -\end{equation} -and apply our initial condition for position, $y(0) = 0$, to obtain $y(t)$.
-Here is the full process: -\begin{align} -\frac{dV}{dt} &= g - \frac{c}{m} V^2 \\ -\frac{dV}{g - \frac{c}{m} V^2} &= dt \\ -\frac{m}{c} \int \frac{dV}{a^2 - V^2} &= \int dt = t + \bar{c}, \quad \text{where} \quad a = \sqrt{\frac{mg}{c}} \\ -\frac{m}{c} \frac{1}{a} \tanh^{-1} \left(\frac{V}{a}\right) &= t + c_1 \\ -V &= a \tanh \left( \frac{a c}{m} t + c_1 \right) \\ -\therefore V(t) &= \sqrt{\frac{mg}{c}} \tanh \left(\sqrt{\frac{gc}{m}} t + c_1\right) -\end{align} -Applying the initial condition for velocity, $V(0) = 0$: -\begin{align} -V(0) &= 0 = \sqrt{\frac{mg}{c}} \tanh \left(0 + c_1\right) \\ -\therefore c_1 &= 0 \\ -V(t) &= \sqrt{\frac{mg}{c}} \tanh \left(\sqrt{\frac{gc}{m}} t\right) -\end{align} -Then, to get $y(t)$, we just need to integrate once more: -\begin{align} -\frac{dy}{dt} = V(t) &= \sqrt{\frac{mg}{c}} \tanh \left(\sqrt{\frac{gc}{m}} t\right) \\ -\int dy &= \sqrt{\frac{mg}{c}} \int \tanh \left(\sqrt{\frac{gc}{m}} t\right) dt \\ -y(t) &= \sqrt{\frac{mg}{c}} \sqrt{\frac{m}{gc}} \log\left(\cosh\left(\sqrt{\frac{gc}{m}} t\right)\right) + c_2 \\ -\rightarrow y(t) &= \frac{m}{c} \log\left(\cosh\left(\sqrt{\frac{gc}{m}} t\right)\right) + c_2 -\end{align} -Finally, we can apply the initial condition for position, $y(0) = 0$, to get our solution: -\begin{align} -y(0) &= 0 = \frac{m}{c} \log\left(\cosh\left(0\right)\right) + c_2 = c_2 \\ -\rightarrow c_2 &= 0 \\ -y(t) &= \frac{m}{c} \log\left(\cosh\left(\sqrt{\frac{gc}{m}} t\right)\right) -\end{align}
- -The catenary problem describes the shape of a hanging chain or rope fixed between two points. (It was also a favorite of one of my professors, Joe Prahl, and I like to teach it as an example in his honor.) The downward displacement of the hanging string/chain/rope as a function of horizontal position, $y(x)$, is governed by the equation -\begin{equation} -y^{\prime\prime} = \sqrt{1 + (y^{\prime})^2} -\end{equation}
-
This is actually a boundary value problem, with the boundary conditions for the displacement at one side $y(0) = 0$ and that the slope is zero in the middle: $\frac{dy}{dx}\left(\frac{L}{2}\right) = 0$. (Please note that I have skipped the derivation of the governing equation, and left some details out.)
-We can solve this via substitution, by letting a new variable $u = y^{\prime}$; then, $u^{\prime} = \frac{du}{dx} = y^{\prime\prime}$. This gives is a first-order ODE, which we can integrate: -\begin{align} -\frac{du}{dx} &= \sqrt{1 + u^2} \\ -\int \frac{du}{\sqrt{1 + u^2}} &= \int dx \\ -\sinh^{-1}(u) &= x + c_1, \quad \text{where } \sinh(x) = \frac{e^x - e^{-x}}{2} \\ -u(x) &= \sinh(x + c_1) -\end{align}
-Then, we can integrate once again to get $y(x)$: -\begin{align} -\frac{dy}{dx} &= u(x) = \sinh(x + c_1) \\ -\int dy &= \int \sinh(x + c_1) dx = \int \left(\sinh(x)\cosh(c_1) + \cosh(x)\sinh(c_1)\right)dx \\ -y(x) &= \cosh(x)\cosh(c_1) + \sinh(x)\sinh(c_1) + c_2 \\ -\rightarrow y(x) &= \cosh(x + c_1) + c_2 -\end{align} -This is the general solution to the catenary problem, and applies to any boundary conditions.
-For our specific case, we can apply the boundary conditions and find the particular solution, though it involves some algebra...: -\begin{align} -y(0) &= 0 = \cosh(c_1) + c_2 \\ -\frac{dy}{dx}\left(\frac{L}{2}\right) &= u(0) = \sinh \left(\frac{L}{2} + c_1\right) \\ -\rightarrow c_1 &= -\frac{L}{2} \\ -0 &= \cosh \left( -\frac{L}{2} \right) + c_2 \\ -\rightarrow c_2 &= -\cosh\left( -\frac{L}{2} \right) = -\cosh\left(\frac{L}{2} \right) -\end{align} -So, the overall solution for the catenary problem with the given boundary conditions is -\begin{equation} -y(x) = \cosh \left(x - \frac{L}{2}\right) - \cosh\left( \frac{L}{2} \right) -\end{equation}
-Let's see what this looks like:
- -L = 1.0;
-x = linspace(0, 1);
-y = cosh(x - (L/2.)) - cosh(L/2.);
-plot(x, y)
-
-Please note that I've made some simplifications in the above work, and skipped the details of how the ODE is derived. In general, the solution for the shape is -\begin{equation} -y(x) = C \cosh \frac{x + c_1}{C} + c_2 -\end{equation} -where you would solve for the constants $C$, $c_1$, and $c_2$ using the constraints: -\begin{align} -\int_{x_a}^{x_b} \sqrt{1 + (y^{\prime})^2} dx &= L \\ -y(x_a) &= y_a \\ -y(x_b) &= y_b \;, -\end{align} -where $L$ is the length of the rope/chain.
-You can read more about the catenary problem here (for example): http://euclid.trentu.ca/aejm/V4N1/Chatterjee.V4N1.pdf
- -An important category of 2nd-order ODEs are those that look like -\begin{equation} -y^{\prime\prime} + p(x) y^{\prime} + q(x) y = 0 -\end{equation} -"Homogeneous" means that the ODE is unforced; that is, the right-hand side is zero.
-Depending on what $p(x)$ and $q(x)$ look like, we have a few different solution approaches:
-First, let's talk about the characteristics of linear, homogeneous 2nd-order ODEs:
-Solutions have two parts: $y(x) = c_1 y_1 + c_2 y_2$, where $y_1$ and $y_2$ are each a basis of the solution.
-The two parts of the solution $y_1$ and $y_2$ are linearly independent.
-One way of defining this is that $a_1 y_1 + a_2 y_2 = 0$ only has the trivial solution $a_1=0$ and $a_2=0$.
-Another way of thinking about this is that $y_1$ and $y_2$ are linearly dependent if one is a multiple of the other, like $y_1 = x$ and $y_2 = 5x$. This cannot be solutions to a linear, homogeneous ODE.
-$y_1$ and $y_2$ each satisfy the ODE. Meaning, you can plug each of them into the ODE for $y$ and obtain 0.
-However, we need both parts together to fully solve the ODE.
-If $y_1$ is known, we can get $y_2$ by reduction of order. Let $y_2 = u y_1$, where $u$ is some unknown function of $x$. Then, put $y_2$ into the ODE $y^{\prime}{\prime} + p(x) y^{\prime} + q(x) y = 0$: -\begin{align} -y_2 &= u y_1 \\ -y_2^{\prime} &= u y_1^{\prime} + u^{\prime} y_1 \\ -y_2^{\prime\prime} &= 2 u^{\prime} y_1^{\prime} + u^{\prime\prime} y_1 + u y_1^{\prime\prime} \\ -\rightarrow u^{\prime\prime} &= - \left[ p(x) + \left(\frac{2 y_1^{\prime}}{y_1}\right) \right] u^{\prime} -\text{or, } u^{\prime\prime} &= - \left( g(x) \right) u^{\prime} -\end{align} -Now, we have an ODE with only $u^{\prime\prime}$, $u^{\prime}$, and some function $g(x)$—so we can solve by substitution! Let $u^{\prime} = v$, and then we have $v^{\prime} = -g(x) v$: -\begin{align} -\frac{dv}{dx} &= - \left( p(x) + \frac{2 y_1^{\prime}}{y_1} \right) v \\ -\int \frac{dv}{v} &= - \int \left(p(x) + \frac{2 y_1^{\prime}}{y_1} \right) dx \\ -\text{Recall } 2 \frac{d}{dx} \left( \ln y_1 \right) &= 2 \frac{y_1^{\prime}}{y_1} \\ -\therefore \int \frac{dv}{v} &= - \int \left(p(x) + 2 \frac{d}{dx} \left( \ln y_1 \right) \right) dx \\ -\ln v &= -\int p(x) dx - 2 \ln y_1 \\ -\rightarrow v &= \frac{\exp\left( -\int p(x)dx \right)}{y_1^2} -\end{align}
-So, the actual solution procedure is then:
-Here's an example, where we know one part of the solution $y_1 = e^{-x}$: -\begin{align} -y^{\prime\prime} + 2 y^{\prime} + y &= 0 \\ -\text{Step 1:} \quad v = \frac{\exp \left( -\int 2dx \right)}{ \left(e^{-x}\right)^2} = \frac{e^{-2x}}{e^{-2x}} &= 1 \\ -\text{Step 2:} \quad u = \int v dx = \int 1 dx &= x \\ -\text{Step 3:} \quad y_2 &= x e^{-x} -\end{align} -Then, the general solution to the ODE is $y(x) = c_1 e^{-x} + c_2 x e^{-x}$.
- -Fourier series are a method we can use to solve inhomogeneous 2nd-order ODEs of the form -\begin{equation} -y^{\prime\prime} + p(t) y^{\prime} + q(t) y = r(t) \;, -\end{equation} -where the forcing function $r(t)$ is periodic. This means looking like one of these examples:
- -% periodic square wave
-subplot(3,1,1);
-squareWave = repmat([1,1,1,1,0,0,0,0], [1, 10]);
-t = linspace(0, 1, length(squareWave));
-plot(t, squareWave); ylim([-0.5, 1.5]);
-
-% sin wave
-subplot(3,1,2);
-t = linspace(0, 1, 100);
-plot(t, sin(t*10*pi)); ylim([-1.5, 1.5]);
-
-% sawtooth
-subplot(3,1,3);
-t = linspace(0, 1, 100);
-y = ((mod(t,2*pi/40)/(pi*2/40))*2)-1;
-plot(t, y); ylim([-1.5, 1.5]);
-
-(Actually, as we'll see later, we can use a Fourier series to represent generic forcing function!)
-Fourier series have been around a while, ever since in 1790 Jean-Baptiste Joseph Fourier found that generic periodic functions could be represented by a sum of series of sin() and cos() functions, harmonically related by frequency.
In general, a Fourier series represents a function $f(t)$ by -\begin{equation} -f(t) = a_0 + \sum_{n=1}^{\infty} a_n \cos (n \omega t) + \sum_{n=1}^{\infty} b_n \sin (n \omega t) -\end{equation} -where $a_0$, $a_n$, and $b_n$ are the Fourier coefficients, $\omega = \frac{2\pi}{T}$ is the frequency of the function $f(t)$, and $T$ is the period. $n$ is an integer used as an index.
- -Considering that $n$ is an integer, and the sine and cosine components of a Fourier series share the same fundamental frequency $\omega$, Fourier series have some useful properties:
-The integral of each component trigonometric function over the period is zero: -\begin{equation} -\int_0^T \sin (n \omega t) dt = 0 = \int_0^T cos (n \omega t) dt -\end{equation}
-The component trigonometric functions are orthogonal over their period: -\begin{equation} -\int_0^T \cos(n \omega t) \sin (m \omega t) dt = 0 -\end{equation} -for all values of $n, m = 1, 2, \ldots, \infty$.
-The component trigonometric functions are also orthogonal with themselves over their period: -\begin{align} -\int_0^T \cos (n \omega t) \cos (m \omega t) dt &= \begin{cases}0 \quad m \neq n \\ \frac{T}{2} \quad m = n \end{cases} \\ -% -\int_0^T \sin (n \omega t) \sin (m \omega t) dt &= \begin{cases}0 \quad m \neq n \\ \frac{T}{2} \quad m = n \end{cases} -\end{align} -for all values of $n, m = 1, 2, \ldots, \infty$.
-We can use the above properties to calculate the Fourier coefficients, given a periodic function $f(t)$. -First, recall -\begin{equation} -f(t) = a_0 + \sum_{n=1}^{\infty} a_n \cos (n \omega t) + \sum_{n=1}^{\infty} b_n \sin (n \omega t) -\end{equation}
-To calculate $a_0$, integrate both sides of the equation over the period: -\begin{equation} -\int_0^T f(t) dt = \int_0^T a_0 dt + \int_0^T \left( \sum_{n=1}^{\infty} a_n \cos (n \omega t) \right)dt + \int_0^T \left( \sum_{n=1}^{\infty} b_n \sin (n \omega t) \right) dt -\end{equation} -For the integrals of the infinite sums, recall that the integral of the sum of some functions is the same as the sum of the integrals of the functions: $\int (a + b + c) = \int a + \int b + \int c$. That means that -\begin{align} -\int_0^T \left( \sum_{n=1}^{\infty} a_n \cos (n \omega t) \right)dt = \int_0^T a_1 \cos (\omega t) dt + \int_0^T a_2 \cos (2 \omega t) dt + \ldots &= 0 \;, \text{ and} \\ -\int_0^T \left( \sum_{n=1}^{\infty} b_n \sin (n \omega t) \right)dt = \int_0^T b_1 \sin (\omega t) dt + \int_0^T b_2 \sin (2 \omega t) dt + \ldots &= 0 \;, -\end{align} -since the integrals of the trigonometric functions are all zero over the period. Thus, -\begin{equation} -a_0 = \frac{1}{T} \int_0^T f(t) dt -\end{equation}
-To calculate $a_n$, multiply both sides of the equation by $\cos(m \omega t)$ and integrate over the period: -\begin{equation} -\int_0^T f(t) \cos(m \omega t) dt = a_0 \int_0^T \cos(m \omega t) dt + \int_0^T \left( \sum_{n=1}^{\infty} a_n \cos (n \omega t) \cos(m \omega t) \right)dt + \int_0^T \left( \sum_{n=1}^{\infty} b_n \sin (n \omega t) \cos(m \omega t) \right) dt -\end{equation} -Let's take a look at each of the three integrals on the right-hand side. First, -\begin{equation} -a_0 \int_0^T \cos(m \omega t) dt = 0 -\end{equation} -because it just integrates cosine over the period. -Skipping to the last term, -\begin{equation} -\int_0^T \left( \sum_{n=1}^{\infty} b_n \sin (n \omega t) \cos(m \omega t) \right) dt = b_1 \int_0^T sin(\omega t) \cos (m \omega t) dt + b_2 \int_0^T \sin (2 \omega t) \cos (m \omega t) dt + \ldots = 0 -\end{equation} -due to orthogonality. We are just left with the middle integral; let's expand a few terms to see what that looks like: -\begin{equation} -\int_0^T \left( \sum_{n=1}^{\infty} a_n \cos (n \omega t) \cos(m \omega t) \right)dt = a_1 \int_0^T \cos (\omega t) \cos (m \omega t) dt + a_2 \int_0^T \cos (2 \omega t) \cos (m \omega t) dt + \ldots -\end{equation} -Again, due to orthogonality, all of the terms of this infinite sum of integrals will be zero, except for the term where $n = m$. As a result, we are left with -\begin{align} -\int_0^T f(t) \cos(m \omega t) dt &= \int_0^T a_m \cos^2 (m \omega t) dt = a_m \frac{T}{2} \\ -a_n = a_m &= \frac{2}{T} \int_0^T f(t) \cos (n \omega t) dt -\end{align}
-We can find $b_n$ in the same way, multiplying the equation by $\sin (m \omega t)$ and integrating over the period. The work is the same, so let's skip that: -\begin{equation} -b_n = \frac{2}{T} \int_0^T f(t) \sin (n \omega t) dt -\end{equation}
-x = [0 1 1 2 2 3 3 4]; y = [2 2 1 1 2 2 1 1];
-plot(x,y); ylim([0 2.5]); ylabel('f(t)'); xlabel('t')
-
-First, we need to identify the fundamental period and frequency: $T = 2$ and then $\omega = \frac{2\pi}{T} = \pi$. Our work is then to calculate the Fourier coefficients. Since our periodic function $f(t)$ is not easily expressed with a function–hence the need for a Fourier series—we'll use piecewise integration.
-First, calculate $a_0$: -\begin{align} -a_0 =& \frac{1}{T} \int_0^T f(t) dt \\ -&= \frac{1}{2} \int_0^2 f(t) dt = \frac{1}{2}\left( \int_0^1 2dt + \int_1^2 1dt \right) = \frac{1}{2} (2\times 1 + 1 \times 1) \\ -a_0 &= \frac{3}{2} -\end{align}
-Then, get $a_n$: -\begin{align} -a_n &= \frac{2}{T} \int_0^T f(t) \cos (n \omega t) dt \\ -&= \frac{2}{2} \int_0^2 f(t) \cos (n \pi t) dt = \left( \int_0^1 2 \cos (n \pi t) dt + \int_1^2 1 \cos (n \pi t)dt \right) \\ -&= \frac{2}{n \pi} \sin(n \pi t)\Big|_0^1 + \frac{1}{n\pi} \sin(n \pi t)\Big|_1^2 \\ -&= \frac{2}{n \pi}\left(\sin(n\pi) - \sin(0)\right) + \frac{1}{n\pi}\left( sin(2n\pi) - \sin(n\pi)\right) \\ -a_n &= 0 -\end{align}
-Finally, we can calculate $b_n$: -\begin{align} -b_n &= \frac{2}{T} \int_0^T f(t) \sin (n \omega t) dt \\ -&= \frac{2}{2} \int_0^2 f(t) \sin (n \pi t) dt = \left( \int_0^1 2 \sin (n \pi t) dt + \int_1^2 1 \sin (n \pi t)dt \right) \\ -&= -\frac{2}{n \pi} \cos(n \pi t)\Big|_0^1 - \frac{1}{n\pi} \cos(n \pi t)\Big|_1^2 \\ -&= -\frac{2}{n \pi}\left(\cos(n\pi) - \cos(0)\right) - \frac{1}{n\pi}\left( cos(2n\pi) - \cos(n\pi)\right) \\ -b_n &= -\frac{2}{n \pi}\left(\cos(n\pi) - 1\right) - \frac{1}{n\pi}\left( 1 - \cos(n\pi)\right) = -\frac{1}{n\pi}\left( \cos(n\pi) - 1\right) -\end{align} -but recall that $n = 1, 2, \ldots, \infty$. Then, -\begin{align} -b_n &= -\frac{1}{n\pi} \times \begin{cases} -2 \text{ if } n \text{ odd} \\0 \text{ if } n \text{ even}\end{cases} \\ -\rightarrow b_n &= \frac{2}{n\pi} \quad n = \text{odd} -\end{align}
-Then, our Fourier series representation for the function shown above is -\begin{equation} -f(t) = \frac{3}{2} + \sum_{\substack{n = 1 \\n = \text{odd}}}^{\infty} \frac{2}{n\pi} \sin (n \pi t) -\end{equation}
-Now, let's see how whether this actually works! Let's start with one term of the infinite sum, then gradually increase.
- -t = linspace(0, 4);
-
-% maximum number of terms
-n_max = [1, 2, 3, 5, 10, 25, 50, 250, 500];
-
-for i = 1 : length(n_max)
- N = n_max(i);
-
- s = 3./2.;
- for n = 1 : 2 : 2*N
- s = s + (2. / (n*pi)) .* sin(n * pi * t);
- end
- subplot(3, 3, i)
- plot(t, s); axis([0 4 0 3]); title(sprintf('%d terms', N));
-end
-
-As we increase the number of terms, adding higher-frequency sine waves, we are better able to match the original rectangular wave. Notice the discrepancies that remain near the sharp corners even after the rest of the series closely resembles the function: these are known as Gibbs phenomena, caused by the Fourier series overshooting or undershooting (or "ringing") near discontinuities.
- -y = [-1 0 1 0 -1];
-x = [-2 -1 0 1 2];
-subplot(2,1,1); plot(x,y); title('Even function');
-ax = gca; ax.XAxisLocation = 'origin'; ax.YAxisLocation = 'origin';
-
-y = [0 -1 0 1 0];
-subplot(2,1,2); plot(x,y); title('Odd function');
-ax = gca; ax.XAxisLocation = 'origin'; ax.YAxisLocation = 'origin';
-
-An even function's Fourier series simplifies to a Fourier cosine series: -\begin{align} -f(x) &= a_0 + \sum_{n=1}^{\infty} a_n \cos (n \omega x) dx \\ -a_0 &= \frac{2}{T} \int_0^{T/2} f(x) dx \\ -a_n &= \frac{4}{T} \int_0^{T/2} f(x) \cos(n \omega x) dx -\end{align}
-An odd function's Fourier series simplifies to a Fourier sine series: -\begin{align} -f(x) &= \sum_{n=1}^{\infty} b_n \sin (n \omega x) dx \\ -b_n &= \frac{4}{T} \int_0^{T/2} f(x) \sin(n \omega x) dx -\end{align}
-Note: not all periodic functions can be considered an even or an odd function.
- -One way we might use a Fourier series is to solve an inhomogeneous 2nd-order ODE, where the forcing term is given by a periodic function not easily expressed using our typical functions.
-For example, let's consider an undamped mass-spring system, where the forcing is given by a periodic rectangular wave: -\begin{equation} -y^{\prime\prime} + 4y = f(t) -\end{equation} -where the forcing function $f(t)$ is
- -t = [0 0 1 1 2 2 3 3];
-y = [0 1 1 -1 -1 1 1 0];
-plot(t,y); ylim([-1.5 1.5]);
-title('Periodic rectangular forcing function');
-ax = gca; ax.XAxisLocation = 'origin';
-
-Using recognizing this as an odd function, we could construct a Fourier sine series to represent the forcing function: -\begin{equation} -f(t) = \sum_{\substack{n=1\\n=\text{odd}}}^{\infty} \frac{4}{n\pi} \sin(n \pi t) -\end{equation}
-Let's confirm this works:
- -t = linspace(0, 3, 500);
-s = 0;
-for n = 1 : 2 : 1000
- s = s + (4/(n*pi)).*sin(n*pi*t);
-end
-plot(t, s)
-
-Looks good!
-To find the exact solution for our displacement $y(t)$, we can follow our usual analytical solution approach: find the homogeneous solution $y_H$, then find the inhomogeneous solution $y_{IH}$; the overall solution is then $y(t) = y_H + y_{IH}$. The homogeneous solution is -\begin{equation} -y_H = c_1 \sin (2t) + c_2 \cos (2t) -\end{equation} -We then find the inhomogeneous solution using -\begin{equation} -y^{\prime\prime} + 4y = \frac{4}{n\pi} \sin (n \pi t) \quad n = 1, 3, \ldots, \infty -\end{equation} -Solving this will give us a specific $y_{IH, n}$; the complete inhomogeneous solution is then -\begin{equation} -y_{IH} = \sum_{\substack{n=1\\n=\text{odd}}}^{\infty} y_{IH, n} \;. -\end{equation}
-Recognizing that our forcing function is sinusoidal, we should use the method of undetermined coefficients: -\begin{equation} -y_{IH, n} = K_1 \sin (n \pi t) + K_2 \cos (n \pi t) -\end{equation} -Inserting this into the above ODE gives -\begin{align} -K_1 &= \frac{4}{n \pi (4 - n^2 \pi^2)} \\ -K_2 &= 0 -\end{align}
-Thus, the overall solution is -\begin{equation} -y(t) = c_1 \sin(2t) + c_2 \cos(2t) + \sum_{\substack{n=1\\n=\text{odd}}}^{\infty} \frac{4}{n \pi (4 - n^2 \pi^2)} \sin (n \pi t) -\end{equation}
- -What about a damped mass-spring system? Recall that the homogeneous solution could take one of these three forms: -\begin{align} -y_H &= c_1 e^{-\lambda_1 t} + c_2 e^{-\lambda_2 t} \\ -y_H &= c_1 e^{-\lambda_1 t} + c_2 t e^{-\lambda_2 t} \text{ or} \\ -y_H &= e^{-\alpha t} (c_1 \sin(\beta t) + c_2 \cos(\beta t)) -\end{align} -while the inhomogeneous solution, given a Fourier series forcing function, will take the form -\begin{equation} -y_{IH} = K_1 \sin() + K_2 \cos() -\end{equation}
-The overall solution combines the homogenenous and inhomogeneous solutions. But, the homogeneous solution in this case is transient, because it decays to zero. On the other hand, the inhomogeneous solution remains, and is the steady-state solution.
- -This section focuses on analytical solutions for initial-value problems, meaning problems where we know the values of $y(t)$ and $\frac{dy}{dt}$ at $t=0$ (or $y(x)$ and $\frac{dy}{dx}$ at $x=0$): $y(0)$ and $y^{\prime}(0)$.
- -An important category of 2nd-order ODEs are those that look like -\begin{equation} -y^{\prime\prime} + p(x) y^{\prime} + q(x) y = 0 -\end{equation} -"Homogeneous" means that the ODE is unforced; that is, the right-hand side is zero.
-Depending on what $p(x)$ and $q(x)$ look like, we have a few different solution approaches:
-First, let's talk about the characteristics of linear, homogeneous 2nd-order ODEs:
-Solutions have two parts: $y(x) = c_1 y_1 + c_2 y_2$, where $y_1$ and $y_2$ are each a basis of the solution.
-The two parts of the solution $y_1$ and $y_2$ are linearly independent.
-One way of defining this is that $a_1 y_1 + a_2 y_2 = 0$ only has the trivial solution $a_1=0$ and $a_2=0$.
-Another way of thinking about this is that $y_1$ and $y_2$ are linearly dependent if one is a multiple of the other, like $y_1 = x$ and $y_2 = 5x$. This cannot be solutions to a linear, homogeneous ODE.
-$y_1$ and $y_2$ each satisfy the ODE. Meaning, you can plug each of them into the ODE for $y$ and obtain 0.
-However, we need both parts together to fully solve the ODE.
-If $y_1$ is known, we can get $y_2$ by reduction of order. Let $y_2 = u y_1$, where $u$ is some unknown function of $x$. Then, put $y_2$ into the ODE $y^{\prime}{\prime} + p(x) y^{\prime} + q(x) y = 0$: -\begin{align} -y_2 &= u y_1 \\ -y_2^{\prime} &= u y_1^{\prime} + u^{\prime} y_1 \\ -y_2^{\prime\prime} &= 2 u^{\prime} y_1^{\prime} + u^{\prime\prime} y_1 + u y_1^{\prime\prime} \\ -\rightarrow u^{\prime\prime} &= - \left[ p(x) + \left(\frac{2 y_1^{\prime}}{y_1}\right) \right] u^{\prime} -\text{or, } u^{\prime\prime} &= - \left( g(x) \right) u^{\prime} -\end{align} -Now, we have an ODE with only $u^{\prime\prime}$, $u^{\prime}$, and some function $g(x)$—so we can solve by substitution! Let $u^{\prime} = v$, and then we have $v^{\prime} = -g(x) v$: -\begin{align} -\frac{dv}{dx} &= - \left( p(x) + \frac{2 y_1^{\prime}}{y_1} \right) v \\ -\int \frac{dv}{v} &= - \int \left(p(x) + \frac{2 y_1^{\prime}}{y_1} \right) dx \\ -\text{Recall } 2 \frac{d}{dx} \left( \ln y_1 \right) &= 2 \frac{y_1^{\prime}}{y_1} \\ -\therefore \int \frac{dv}{v} &= - \int \left(p(x) + 2 \frac{d}{dx} \left( \ln y_1 \right) \right) dx \\ -\ln v &= -\int p(x) dx - 2 \ln y_1 \\ -\rightarrow v &= \frac{\exp\left( -\int p(x)dx \right)}{y_1^2} -\end{align}
-So, the actual solution procedure is then:
-Here's an example, where we know one part of the solution $y_1 = e^{-x}$: -\begin{align} -y^{\prime\prime} + 2 y^{\prime} + y &= 0 \\ -\text{Step 1:} \quad v = \frac{\exp \left( -\int 2dx \right)}{ \left(e^{-x}\right)^2} = \frac{e^{-2x}}{e^{-2x}} &= 1 \\ -\text{Step 2:} \quad u = \int v dx = \int 1 dx &= x \\ -\text{Step 3:} \quad y_2 &= x e^{-x} -\end{align} -Then, the general solution to the ODE is $y(x) = c_1 e^{-x} + c_2 x e^{-x}$.
- -A common category of 2nd-order homogeneous ODEs are equations with constant coefficients, of the form: -\begin{equation} -y^{\prime\prime} + a y^{\prime} + by = 0 -\end{equation} -Note that these are unforced, and the right-hand side is zero.
-Solutions to these equations take the form $y(x) = e^{\lambda x}$, and inserting this into the ODE gives us the characteristic equation -\begin{equation} -\lambda^2 + a \lambda + b = 0 -\end{equation} -which we can solve to find the solution for given coefficients $a$ and $b$ and initial conditions. Depending on those coefficients and the solution to the characteristic equation, our solution can fall into one of three cases:
-Real roots: $\lambda_1$ and $\lambda_2$. This is an overdamped system and the full solution takes the form -\begin{equation} -y(x) = c_1 e^{\lambda_1 x} + c_2 e^{\lambda_2 x} -\end{equation}
-Repeated roots: $\lambda_1 = \lambda_2 = \lambda$. This is a critically damped system and the full solution is -\begin{equation} -y(x) = c_1 e^{\lambda x} + c_2 x e^{\lambda x} -\end{equation} -(Where does that second part come from, you might ask? Well, we know that $y_1$ is $e^{\lambda x}$, but the second part cannot also be $e^{\lambda x}$ because those are linearly dependent. So, we use reduction of order to find $y_2$, which is $x e^{\lambda x}$.
-Imaginary roots: $\lambda = \frac{-a}{2} \pm \beta i$, where $\beta = \frac{1}{2} \sqrt{4b - a^2}$. This is an underdamped system and the solution takes the form -\begin{equation} -y(x) = e^{-ax/2} \left( c_1 \sin \beta x + c_2 \cos \beta x \right) -\end{equation}
-Some examples:
-$y^{\prime\prime} + 3 y^{\prime} + 2y = 0$ -\begin{align} -\rightarrow \lambda^2 + 3\lambda + 2 &= 0 \\ -(\lambda + 2)(\lambda + 1) &= 0 \\ -\lambda &= -2, -1 \\ -y(x) &= c_1 e^{-x} + c_2 e^{-2x} -\end{align} -Then, we would use the initial conditions given for $y(0)$ and $y^{\prime}(0)$ to find $c_1$ and $c_2$.
-$y^{\prime\prime} + 6 y^{\prime} + 9y = 0$ -\begin{align} -\rightarrow \lambda^2 + 6\lambda + 9 &= 0 \\ -(\lambda + 3)(\lambda + 3) &= 0 \\ -\lambda &= -3 \\ -y(x) &= c_1 e^{-3x} + c_2 x e^{-3x} -\end{align}
-$y^{\prime\prime} + 6 y^{\prime} + 25 y = 0$ -\begin{align} -\rightarrow \lambda^2 + 6\lambda + 25 &= 0 \\ -\lambda &= -3 \pm 4i \\ -y(x) &= e^{-3x} \left( c_1 \sin 4x + c_2 \cos 4x \right) -\end{align}
-Euler-Cauchy equations are of the form -\begin{equation} -x^2 y^{\prime\prime} + axy^{\prime} + by = 0 -\end{equation}
-Solutions take the form $y = x^m$, which when plugged into the ODE leads to a different characterisic equation to find $m$: -\begin{align} -y &= x^m \\ -y^{\prime} &= m x^{m-1} \\ -y^{\prime\prime} &= m (m-1) x^{m-2} \\ -\rightarrow x^2 m (m-1) x^{m-2} + axmx^{m-1} + bx^m &= 0 \\ -m^2 + (a-1)m + b &= 0 -\end{align} -This is our new characteristic formula for these problems, and solving for the roots of this equation gives us $m$ and thus our general solution.
-Like equations with constant coefficients, we have three solution forms depending on the roots of the characteristic equation:
-Inhomogeneous, or forced, 2nd-order ODEs with constant coefficients take the form -\begin{equation} -y^{\prime\prime} + a y^{\prime} + by = F(t) -\end{equation} -with initial conditions $y(0) = y_0$ and $y^{\prime}(0) = y_0^{\prime}$. Depending on the form of the forcing function $F(t)$, we can solve with techniques such as
-The solution in general to inhomogeneous ODEs includes two parts: -\begin{equation} -y(t) = y_{\text{H}} + y_{\text{IH}} = c_1 y_1 + c_2 y_2 + y_{\text{IH}} \;, -\end{equation} -where $y_{\text{H}}$ is the solution from the equivalent homogeneous ODE $y^{\prime\prime} + a y^{\prime} + b y = 0$.
-The forcing function $F(t)$ may be
-For continuous forcing functions, we have two solution methods: the method of undetermined coefficients, and variation of parameters.
-Generally you'll want to use the method of undetermined coefficients when possible, which depends on if $F(t)$ matches one of a set of functions. In that case, the form of the inhomogeneous solution $y_{\text{IH}}(t)$ follows that of the forcing function $F(t)$, with one or more unknown constants:
-| $F(t)$ | -$y_{\text{IH}}(t)$ | -
|---|---|
| constant | -$K$ | -
| $\cos \omega t$ | -$K_1 \cos \omega t + K_2 \sin \omega t$ | -
| $\sin \omega t$ | -$K_1 \cos \omega t + K_2 \sin \omega t$ | -
| $e^{-at}$ | -$K e^{-at}$ | -
| $(A) t$ | -$K_0 + K_1 t$ | -
| $t^n$ | -$K_0 + K_1 t + K_2 t^2 + \ldots + K_n t^n$ | -
For combinations of these functions, we can combine functions; for example, given -\begin{align} -F(t) &= e^{-at} \cos \omega t \quad \text{or} e^{-at} \sin \omega t \\ -y_{\text{IH}} &= K_1 e^{-at} \cos \omega t + K_2 e^{-at} \sin \omega t -\end{align} -(Note how in all the above cases how the inhomogeneous solution follows the functional form of the forcing function; for example, the exponential decay rate $a$ or the sinusoidal frequency $\omega$ match.
-The method of undetermined coefficients works by plugging the candidate inhomogeneous solutionn $y_{\text{IH}}$ into the full ODE, and solving for the constants (e.g., $K$)—but not from the initial conditions.
-For example, let's solve -\begin{equation} -y^{\prime\prime} + 2y^{\prime} + y = e^{-x} -\end{equation} -with initial conditions $y(0) = y^{\prime}(0) = 0$. First, we should find the solution to the homogeneous equation -\begin{equation} -y^{\prime\prime} + 2y^{\prime} + y = 0 \;. -\end{equation} -We can do this by using the associated characteristic formula -\begin{align} -\lambda^2 + 2 \lambda + 1 &= 0 \\ -(\lambda + 1)(\lambda + 1) &= 0 \\ -\rightarrow y_{\text{H}} &= c_1 e^{-x} + c_2 x e^{-x} -\end{align}
-To find the inhomogeneous solution, we would look at the table above to find what matches the forcing function $e^{-x}$. Normally, we'd grab $K e^{-x}$, but that would not be linearly independent from the first part of the homogeneous solution $y_{\text{H}}$. The same is true for $K x e^{-x}$, which is linearly dependent with the second part of $y_{\text{H}}$, but $K x^2 e^{-x}$ works! Then, we just need to find $K$ by plugging this into the ODE: -\begin{align} -y_{\text{IH}} &= K x^2 e^{-x} \\ -y^{\prime} &= K e^{-x} (2x - x^2) \\ -y^{\prime\prime} &= K e^{-x} (x^2 - 4x + 2) \\ -2 K &= 1 \\ -\rightarrow K &= \frac{1}{2} \\ -y_{\text{IH}} &= \frac{1}{2} x^2 e^{-x} -\end{align} -Thus, the overall general solution is -\begin{equation} -y(x) = c_1 e^{-x} + c_2 x e^{-x} + \frac{1}{2} x^2 e^{-x} -\end{equation} -and we would solve for the integration constants $c_1$ and $c_2$ using the initial conditions.
-Important points to remember:
-We have the variation of parameters approach to solve for inhomogeneous 2nd-order ODEs that are more general: -\begin{equation} -y^{\prime\prime} + p(x) y^{\prime} + q(x) y = r(x) -\end{equation} -In this case, we can assume a solution $y(x) = y_1 u_1 + y_2 u_2$.
-The solution procedure is:
-Obtain $y_1$ and $y_2$ by solving the homogeneous equation: $y^{\prime\prime} + p(x) y^{\prime} + q(x) y = 0$
-Solve for $u_1$ and $u_2$: -\begin{align} -u_1 &= - \int \frac{y_2 r(x)}{W} dx + c_1 \\ -u_2 &= \int \frac{y_1 r(x)}{W} dx + c_2 \\ -W &= \begin{vmatrix} -y_1 & y_2\\ y_1^{\prime} & y_2^{\prime}\\ -\end{vmatrix} = y_1 y_2^{\prime} - y_2 y_1^{\prime} \;, -\end{align} -where $W$ is the Wronksian.
-Then, we have the general solution: -\begin{align} -y &= u_1 y_1 + u_2 y_2 \\ -&= \left( -\int \frac{y_2 r(x)}{W} dx + c_1 \right) y_1 + \left( \int \frac{y_1 r(x)}{W} dx + c_2 \right) y_2 \;, -\end{align} -where we solve for $c_1$ and $c_2$ using the two initial conditions.
-First, let's try the same example we used for the method of undetermined coefficients above: -\begin{equation} -y^{\prime\prime} + 2 y^{\prime} + y = e^{-x} -\end{equation} -We already found the homogeneous solution, so we know that $y_1 = e^{-x}$ and $y_2 = x e^{-x}$. -Next, let's get the Wronksian, and then $u_1$ and $u_2$. -\begin{align} -W &= \begin{vmatrix} y_1 & y_2 \\ y_1^{\prime} & y_2^{\prime} \end{vmatrix} = e^{-x} e^{-x}(1-x) - x e^{-x} (-e^{-x}) = e^{-2x} \\ -% -u_1 &= -\int \frac{x e^{-x} e^{-x}}{e^{-2x}} dx + c_1 = -\int x dx + c_1 = -\frac{1}{2} x^2 + c_1 \\ -u_2 &= \int \frac{e^{-x} e^{-x}}{e^{-2x}} dx + c_2 = \int dx + c_2 = x + c_2 \\ -y(x) &= \left(-\frac{1}{2} x^2 + c_1\right) e^{-x} + (x + c_2) x e^{-x} \\ -\end{align} -After simplifying, we obtain the same solution as via the method of undetermined coefficients (but with a bit more work): -\begin{equation} -y(x) = x_1 e^{-x} + c_2 x e^{-x} + \frac{1}{2} x^2 e^{-x} -\end{equation}
-Now let's try an example that we could not solve using the method of undetermined coefficients, with a forcing term that involves hyperbolic cosine (cosh); recall that $\cosh(x) = \frac{e^x + e^{-x}}{2}$. -\begin{equation} -y^{\prime\prime} + 4 y^{\prime} + 4y = \cosh(x) -\end{equation} -First, we need to find the homogeneous solution: -\begin{align} -y^{\prime\prime} + 4 y^{\prime} + 4y &= 0 \\ -\lambda^2 + 4 \lambda + 4 &= 0 \\ -\rightarrow \lambda &= -2 -\end{align} -So our homogeneous solution involves repeated roots: -\begin{equation} -y_H = c_1 e^{-2x} + c_2 x e^{-2x} -\end{equation} -where $y_1 = e^{-2x}$ and $y_2 = x e^{-2x}$.
-Then, we need to find $u_1$ and $u_2$, so let's get the Wronksian and then solve -\begin{align} -W &= \begin{vmatrix} y_1 & y_2 \\ y_1^{\prime} & y_2^{\prime} \end{vmatrix} = e^{-2x} (e^{-2x}) (1 - 2x) - x e^{-2x}(-2 e^{-2x}) = e^{-4x} \\ -% -u_1 &= - \int \frac{x e^{-2x} \cosh x}{e^{-4x}} dx + c_1 = -\int \frac{x \frac{1}{2}(e^x + e^{-x})}{e^{-2x}} dx + c_1 \\ - &= -\frac{1}{2} \int x (e^{3x} + e^x) dx + c_1 = -\frac{1}{2} \left[ \frac{1}{9} e^{3x}(3x-1) + e^x(x-1) \right] + c_1 \\ -u_1 &= -\frac{1}{18} e^{3x}(3x-1) - \frac{1}{2} e^x (x-1) + c_1 \\ -% -u_2 &= \int \frac{e^{-2x} \cosh x}{e^{-4x}} dx + c_2 = \frac{1}{2} \int e^{2x}(e^x + e^{-x}) dx + c_2 = \frac{1}{2} \int (e^{3x} + e^x) dx + c_2 \\ -u_2 &= \frac{1}{6} e^{3x} + \frac{1}{2} e^x + c_2 -\end{align}
-Then, when we put these all together, we get the full (complicated) solution: -\begin{equation} -y(x) = \left[ -\frac{1}{18} e^{3x} (3x-1) - \frac{1}{2} e^x (x-1) + c_1 \right] e^{-2x} + \left( \frac{1}{6} e^{3x} + \frac{1}{2} e^x + c_2 \right) x e^{-2x} -\end{equation}
- -We've gone over how to solve 1st-order ODEs using numerical methods, but what about 2nd-order or any higher-order ODEs? We can use the same methods we've already discussed by transforming our higher-order ODEs into a system of first-order ODEs.
-Meaning, if we are given a 2nd-order ODE -\begin{equation} -\frac{d^2 y}{dx^2} = y^{\prime\prime} = f(x, y, y^{\prime}) -\end{equation} -we can transform this into a system of two 1st-order ODEs that are coupled: -\begin{align} -\frac{dy}{dx} &= y^{\prime} = u \\ -\frac{du}{dx} &= u^{\prime} = y^{\prime\prime} = f(x, y, u) -\end{align} -where $f(x, y, u)$ is the same as that given above for $\frac{d^2 y}{dx^2}$.
-Thus, instead of a 2nd-order ODE to solve, we have two 1st-order ODEs: -\begin{align} -y^{\prime} &= u \\ -u^{\prime} &= f(x, y, u) -\end{align}
-So, we can use all of the methods we have talked about so far to solve 2nd-order ODEs by transforming the one equation into a system of two 1st-order equations.
-This works for higher-order ODEs too! For example, if we have a 3rd-order ODE, we can transform it into a system of three 1st-order ODEs: -\begin{align} -\frac{d^3 y}{dx^3} &= f(x, y, y^{\prime}, y^{\prime\prime}) \\ -\rightarrow y^{\prime} &= u \\ -u^{\prime} &= y^{\prime\prime} = w \\ -w^{\prime} &= y^{\prime\prime\prime} = f(x, y, u, w) -\end{align}
- -For example, let's solve a forced damped mass-spring problem given by a 2nd-order ODE: -\begin{equation} -y^{\prime\prime} + 5y^{\prime} + 6y = 10 \sin \omega t -\end{equation} -with the initial conditions $y(0) = 0$ and $y^{\prime}(0) = 5$.
-We start by transforming the equation into two 1st-order ODEs. Let's use the variables $z_1 = y$ and $z_2 = y^{\prime}$: -\begin{align} -\frac{dz_1}{dt} &= z_1^{\prime} = z_2 \\ -\frac{dz_2}{dt} &= z_2^{\prime} = y^{\prime\prime} = 10 \sin \omega t - 5z_2 - 6z_1 -\end{align}
-Then, let's solve numerically using the forward Euler method. Recall that the recursion formula for forward Euler is: -\begin{equation} -y_{i+1} = y_i + \Delta x f(x_i, y_i) -\end{equation} -where $f(x,y) = \frac{dy}{dx}$.
-Let's solve using $\omega = 1$ and with a step size of $\Delta t = 0.1$, over $0 \leq t \leq 3$.
-We can compare this against the exact solution, obtainable using the method of undetermined coefficients: -\begin{equation} -y(t) = -6 e^{-3t} + 7 e^{-2t} + \sin t - \cos t -\end{equation}
- -% plot exact solution first
-t = linspace(0, 3);
-y_exact = -6*exp(-3*t) + 7*exp(-2*t) + sin(t) - cos(t);
-plot(t, y_exact); hold on
-
-omega = 1;
-
-dt = 0.1;
-t = [0 : dt : 3];
-
-f = @(t,z1,z2) 10*sin(omega*t) - 5*z2 - 6*z1;
-
-z1 = zeros(length(t), 1);
-z2 = zeros(length(t), 1);
-z1(1) = 0;
-z2(1) = 5;
-for i = 1 : length(t)-1
- z1(i+1) = z1(i) + dt * z2(i);
- z2(i+1) = z2(i) + dt * f(t(i), z1(i), z2(i));
-end
-
-plot(t, z1, 'o--')
-xlabel('time'); ylabel('displacement')
-legend('Exact', 'Forward Euler', 'Location','southeast')
-
-clear
-% plot exact solution first
-t = linspace(0, 3);
-y_exact = -6*exp(-3*t) + 7*exp(-2*t) + sin(t) - cos(t);
-plot(t, y_exact); hold on
-
-omega = 1;
-
-dt = 0.1;
-t = [0 : dt : 3];
-
-f = @(t,z1,z2) 10*sin(omega*t) - 5*z2 - 6*z1;
-
-z1 = zeros(length(t), 1);
-z2 = zeros(length(t), 1);
-z1(1) = 0;
-z2(1) = 5;
-for i = 1 : length(t)-1
- % predictor
- z1p = z1(i) + z2(i)*dt;
- z2p = z2(i) + f(t(i), z1(i), z2(i))*dt;
-
- % corrector
- z1(i+1) = z1(i) + 0.5*dt*(z2(i) + z2p);
- z2(i+1) = z2(i) + 0.5*dt*(f(t(i), z1(i), z2(i)) + f(t(i+1), z1p, z2p));
-end
-plot(t, z1, 'o')
-xlabel('time'); ylabel('displacement')
-legend('Exact', 'Heuns', 'Location','southeast')
-
-ode45 We can also solve using ode45, by providing a separate function file that defines the system of 1st-order ODEs. In this case, we'll need to use a single array variable, Z, to store $z_1$ and $z_2$. The first column of Z will store $z_1$ (Z(:,1)) and the second column will store $z_2$ (Z(:,2)).
%%file mass_spring.m
-function dzdt = mass_spring(t, z)
- omega = 1;
- dzdt = zeros(2,1);
-
- dzdt(1) = z(2);
- dzdt(2) = 10*sin(omega*t) - 6*z(1) - 5*z(2);
-end
-% plot exact solution first
-t = linspace(0, 3);
-y_exact = -6*exp(-3*t) + 7*exp(-2*t) + sin(t) - cos(t);
-plot(t, y_exact); hold on
-
-% solution via ode45:
-[T, Z] = ode45('mass_spring', [0 3], [0 5]);
-
-plot(T, Z(:,1), 'o')
-xlabel('time'); ylabel('displacement')
-legend('Exact', 'ode45', 'Location','southeast')
-
-We saw how to implement the Backward Euler method for a 1st-order ODE, but what about a 2nd-order ODE? (Or in general a system of 1st-order ODEs?)
-The recursion formula is the same, except now our dependent variable is an array/vector: -\begin{equation} -\mathbf{y}_{i+1} = \mathbf{y}_i + \Delta t \, \mathbf{f} \left( t_{i+1} , \mathbf{y}_{i+1} \right) -\end{equation} -where the bolded $\mathbf{y}$ and $\mathbf{f}$ indicate array quantities (in other words, they hold more than one value).
-In general, we can use Backward Euler to solve 2nd-order ODEs in a similar fashion as our other numerical methods:
-The main difference is that we will now have a system of two equations and two unknowns: $y_{1, i+1}$ and $y_{2, i+1}$.
-Let's demonstrate with an example: -\begin{equation} -y^{\prime\prime} + 6 y^{\prime} + 5y = 10 \quad y(0) = 0 \quad y^{\prime}(0) = 5 -\end{equation} -where the exact solution is -\begin{equation} -y(t) = -\frac{3}{4} e^{-5t} - \frac{5}{4} e^{-t} + 2 -\end{equation}
-To solve numerically,
-Convert the 2nd-order ODE into a system of two 1st-order ODEs: -\begin{gather} -y_1 = y \quad y_1(t=0) = 0 \\ -y_2 = y^{\prime} \quad y_2 (t=0) = 5 -\end{gather} -Then, for the derivatives of these variables: -\begin{align} -y_1^{\prime} &= y_2 \\ -y_2^{\prime} &= 10 - 6 y_2 - 5 y_1 -\end{align}
-Then plug these derivatives into the Backward Euler recursion formulas and solve for $y_{1,i+1}$ and $y_{2,i+1}$: -\begin{align} -y_{1, i+1} &= y_{1, i} + \Delta t \, y_{2,i+1} \\ -y_{2, i+1} &= y_{2, i} + \Delta t \left( 10 - 6 y_{2, i+1} - 5 y_{1,i+1} \right) \\ -\\ -y_{1, i+1} - \Delta t \, y_{2, i+1} &= y_{1,i} \\ -5 \Delta t \, y_{1, i+1} + (1 + 6 \Delta t) y_{2, i+1} &= y_{2,i} + 10 \Delta t \\ -\text{or} \quad -\begin{bmatrix} 1 & -\Delta t \\ 5 \Delta t & (1+6\Delta t)\end{bmatrix} -\begin{bmatrix} y_{1, i+1} \\ y_{2, i+1} \end{bmatrix} &= -\begin{bmatrix} y_{1,i} \\ y_{2,i} + 10 \Delta t \\ \end{bmatrix} \\ -\mathbf{A} \mathbf{y}_{i+1} &= \mathbf{b} -\end{align} -To isolate $\mathbf{y}_{i+1}$ and get a usable recursion formula, we need to solve this system of two equations. We could solve this by hand using the substitution method, or we can use Cramer's rule: -\begin{align} -y_{1, i+1} &= \frac{ y_{1,i} (1 + 6 \Delta t) + \Delta t \left( y_{2,i} + 10 \Delta t \right)}{1 + 6 \Delta t + 5 \Delta t^2} \\ -y_{2, i+1} &= \frac{ y_{2,i} + 10 \Delta t - 5 \Delta t y_{1,i}}{1 + 6 \Delta t + 5 \Delta t^2} -\end{align}
-Let's confirm that this gives us a good, well-behaved numerical solution and compare with the Forward Euler method:
- -clear
-
-% Exact solution
-t = linspace(0, 5);
-y_exact = @(t) -(3/4)*exp(-5*t) - (5/4)*exp(-t) + 2;
-plot(t, y_exact(t)); hold on
-
-dt = 0.1;
-t = 0 : dt : 5;
-
-% Forward Euler
-f = @(t, y1, y2) 10 - 6*y2 - 5*y1;
-y1 = zeros(length(t), 1); y2 = zeros(length(t), 1);
-y1(1) = 0; y2(1) = 5;
-for i = 1 : length(t) - 1
- y1(i+1) = y1(i) + dt*y2(i);
- y2(i+1) = y2(i) + dt*f(t(i), y1(i), y2(i));
-end
-plot(t, y1, '+')
-
-Y = zeros(length(t), 2);
-Y(1,1) = 0;
-Y(1,2) = 5;
-for i = 1 : length(t) - 1
- D = 1 + 6*dt + 5*dt^2;
- Y(i+1, 1) = (Y(i,1)*(1 + 6*dt) + dt*(Y(i,2) + 10*dt)) / D;
- Y(i+1, 2) = (Y(i,2) + 10*dt - Y(i,1)*5*dt) / D;
-end
-plot(t, Y(:,1), 'o')
-legend('Exact', 'Forward Euler', 'Backward Euler', 'location', 'southeast')
-
-fprintf('Maximum error of Forward Euler: %5.3f\n', max(abs(y1(:) - y_exact(t)')));
-fprintf('Maximum error of Backward Euler: %5.3f', max(abs(Y(:,1) - y_exact(t)')));
-
-So, for $\Delta t = 0.1$, we see that the Forward and Backward Euler methods give an error $\mathcal{O}(\Delta t)$, as expected since both methods are first-order accurate.
-Let's see how they compare for a larger step size:
- -clear
-
-% Exact solution
-t = linspace(0, 5);
-y_exact = @(t) -(3/4)*exp(-5*t) - (5/4)*exp(-t) + 2;
-plot(t, y_exact(t)); hold on
-
-dt = 0.5;
-t = 0 : dt : 5;
-
-% Forward Euler
-f = @(t, y1, y2) 10 - 6*y2 - 5*y1;
-y1 = zeros(length(t), 1); y2 = zeros(length(t), 1);
-y1(1) = 0; y2(1) = 5;
-for i = 1 : length(t) - 1
- y1(i+1) = y1(i) + dt*y2(i);
- y2(i+1) = y2(i) + dt*f(t(i), y1(i), y2(i));
-end
-plot(t, y1, 'o')
-
-% Backward Euler
-
-Y = zeros(length(t), 2);
-Y(1,1) = 0;
-Y(1,2) = 5;
-for i = 1 : length(t) - 1
- D = 1 + 6*dt + 5*dt^2;
- Y(i+1, 1) = (Y(i,1)*(1 + 6*dt) + dt*(Y(i,2) + 10*dt)) / D;
- Y(i+1, 2) = (Y(i,2) + 10*dt - Y(i,1)*5*dt) / D;
-end
-plot(t, Y(:,1), 'o')
-legend('Exact', 'Backward Euler', 'location', 'southeast')
-
-%fprintf('Maximum error of Forward Euler: %5.3f\n', max(abs(y1(:) - y_exact(t)')));
-%fprintf('Maximum error of Backward Euler: %5.3f', max(abs(Y(:,1) - y_exact(t)')));
-
-fprintf('Maximum error of Forward Euler: %5.3f\n', max(abs(y1(:) - y_exact(t)')));
-fprintf('Maximum error of Backward Euler: %5.3f', max(abs(Y(:,1) - y_exact(t)')));
-
-Backward Euler, since it is unconditionally stable, remains well-behaved at this larger step size, while the Forward Euler method blows up.
-One other thing: instead of using Cramer's rule to get expressions for $y_{1,i+1}$ and $y_{2,i+1}$, we could instead use Matlab to solve the linear system of equations at each time step. To do that, we could replace
-A = [1 -dt; 5*dt (1+6*dt)];
-b = [Y(i,1); Y(i,2)+10*dt];
-Y(i+1,:) = (A\b)';
-where A\b is equivalent to inv(A)*b, but faster. Let's confirm that this gives the same answer:
clear
-
-% Exact solution
-t = linspace(0, 5);
-y_exact = @(t) -(3/4)*exp(-5*t) - (5/4)*exp(-t) + 2;
-plot(t, y_exact(t)); hold on
-
-dt = 0.1;
-t = 0 : dt : 5;
-
-Y = zeros(length(t), 2);
-Y(1,1) = 0;
-Y(1,2) = 5;
-for i = 1 : length(t) - 1
- A = [1 -dt; 5*dt (1+6*dt)];
- b = [Y(i,1); Y(i,2)+10*dt];
- Y(i+1,:) = (A\b)';
-end
-plot(t, Y(:,1), 'o')
-legend('Exact', 'Backward Euler', 'location', 'southeast')
-
-Cramer's Rule provides a solution method for a system of linear equations, where the number of equations equals the number of unknowns. It works for any number of equations/unknowns, but isn't really practical for more than two or three. We'll focus on using it for a system of two equations, with two unknowns $x_1$ and $x_2$: -\begin{gather} -a_{11} + x_1 + a_{12} x_2 = b_1 \\ -a_{21} + x_1 + a_{22} x_2 = b_2 \\ -\text{or } \mathbf{A} \mathbf{x} = \mathbf{b} -\end{gather} -where -\begin{gather} -\mathbf{A} = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \\ -\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\ -\mathbf{b} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix} -\end{gather}
-The solutions for the unknowns are then -\begin{align} -x_1 &= \frac{ \begin{vmatrix} b_1 & a_{12} \\ b_2 & a_{22} \end{vmatrix} }{D} = \frac{b_1 a_{22} - a_{12} b_2}{D} \\ -x_2 &= \frac{ \begin{vmatrix} a_{11} & b_1 \\ a_{21} & b_2 \end{vmatrix} }{D} = \frac{a_{11} b_2 - b_1 a_{21}}{D} -\end{align} -where $D$ is the determinant of $\mathbf{A}$: -\begin{equation} -D = \det(\mathbf{A}) = \begin{vmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{vmatrix} = a_{11} a_{22} - a_{12} a_{21} -\end{equation} -This works as long as the determinant does not equal zero.
- -Power series solutions are another technique we can use to solve 2nd-order homogeneous ODEs of the form -\begin{equation} -y^{\prime\prime} + p(x) y^{\prime} + q(x) y = 0 -\end{equation} -This is useful for more-general cases where our other techniques fail.
-For example, how would you find the solution to this ODE? -\begin{equation} -(1 + x^2) y^{\prime\prime} - 4 x y^{\prime} + 6y = 0 -\end{equation} -None of the methods we've discussed so far would allow us to find an analytical solution to this problem—but we can using a power series solution.
-Power series solutions will be of the form -\begin{equation} -y = \sum_{n=0}^{\infty} a_n x^n -\end{equation} -where the coefficients $a_n$ are what we need to find.
-Continuous means that there should be no discontinuity at $x=0$.
-Bounded means that the solution should be finite at $x=0$.
-For example, consider the ODE -\begin{equation} -y^{\prime\prime} - 4xy^{\prime} + (4x^2 - 2)y = 0 -\end{equation} -Both $p(x) = -4x$ and $q(x) = (4x^2 - 2)$ are continuous and bounded at $x=0$, so $x=0$ is an ordinary point.
-On the other hand, what about -\begin{equation} -y^{\prime\prime} + x^3 y^{\prime} + \frac{1}{x} y = 0 \text{ ?} -\end{equation} -In this case, the solution is unbounded at $x=0$, and so it is not an ordinary point.
-We then solve for the coefficients $a_n$ by plugging this in to the ODE. To do that, we'll need to take advantage of certain properties of power series.
-Dummy index rule. We can replace the index variable used in the power series with another index variable arbitrarily: -\begin{equation} -\sum_{n=0}^{\infty} a_n x^n = \sum_{m=0}^{\infty} a_m x^m -\end{equation} -This is because the index variable is just a "dummy" that only has meaning inside the sum.
-Product rule. We can bring variables, including $x$, multiplying an entire power series into the power series: -\begin{equation} -x \sum_{n=0}^{\infty} a_n x^n = \sum_{n=0}^{\infty} a_n x^{n+1} -\end{equation}
-Derivatives. We can take derivatives of our power series: -\begin{align} -y^{\prime} &= \sum_{n=0}^{\infty} a_n (n) x^{n-1} = \sum_{n=1}^{\infty} a_n (n) x^{n-1} \\ -y^{\prime\prime} &= \sum_{n=0}^{\infty} a_n (n)(n-1) x^{n-2} = \sum_{n=2}^{\infty} a_n (n)(n-1) x^{n-2} -\end{align} -Notice that we can change where the sums in the power series start, because for $y^{\prime}$ the term corresponding to $n=1$ would just be zero, and similar for the first two terms of $y^{\prime\prime}$.
-Index shift. We can redefine the index used within a sum to shift where it starts. For example, if we let $m=n-1$, or $n=m+1$, then: -\begin{equation} -\sum_{n=1}^{\infty} a_n (n) x^{n-1} = \sum_{m=0}^{\infty} a_{m+1} (m+1) x^m -\end{equation} -Or, in other case, if we let $m=n-2$, or $n=m+2$, then: -\begin{equation} -\sum_{n=2}^{\infty} a_n (n)(n-1) x^{n-2} = \sum_{m=0}^{\infty} a_{m+2} (m+2)(m+1) x^m -\end{equation}
-Now, let's apply these properties to solve ODEs.
-Let's try to apply the power series approach to solve -\begin{equation} -y^{\prime\prime} + y = 0 \;, -\end{equation} -where we know the solution will be $y(x) = c_1 \sin x + c_2 \cos x$.
-Is $x=0$ an ordinary point? Yes, the ODE is continuous and bounded at $x=0$. -So, we can find a solution of the form $y(x) = \sum_{n=0}^{\infty} a_n x^n$.
-Now, we solve for the coefficents by plugging the power series into the ODE: -\begin{equation} -\sum_{n=2}^{\infty} a_n (n)(n-1) x^{n-2} + \sum_{n=0}^{\infty} a_n x^n = 0 -\end{equation} -Let's use the index shift rule on the first part of that: -\begin{equation} -\sum_{n=2}^{\infty} a_n (n)(n-1) x^{n-2} \rightarrow \sum_{m=0}^{\infty} a_{m+2} (m+2)(m+1) x^m -\end{equation} -Then, we can use the dummy index rule to change $m$ back to $n$: -\begin{equation} -\sum_{m=0}^{\infty} a_m (m+2)(m+1) x^m \rightarrow \sum_{n=0}^{\infty} a_n (n+2)(n+1) x^n -\end{equation} -Now, let's replace the first term in the ODE with that, merge both terms into a single sum, and simplify: -\begin{align} -\sum_{n=0}^{\infty} a_n (n+2)(n+1) x^n + \sum_{n=0}^{\infty} a_n x^n &= 0 \\ -\sum_{n=0}^{\infty} x^n \left[ a_{n+2}(n+2)(n+1) + a_n \right] &= 0 -\end{align}
-There are infinite terms in this sum, involving the continuous variable $x$; the only way that equation can be satisfied is if
-Use that expression to define a recursive formula for the coefficients: -\begin{equation} -a_{n+2} = \frac{-a_n}{(n+1)(n+2)} -\end{equation} -We can see that the even coefficients will be related to each other, and the odd coefficients will be related. Let's try to identify a pattern with each, starting with the even terms: -\begin{align} -n=0: \quad a_2 &= \frac{-a_0}{1 \cdot 2} = \frac{-a_0}{2!} \\ -n=2: \quad a_4 &= \frac{-a_2}{3 \cdot 4} = \frac{a_0}{4!} \\ -n=4: \quad a_6 &= \frac{-a_4}{5 \cdot 6} = \frac{-a_0}{6!} -\end{align} -and the odd terms: -\begin{align} -n=1: \quad a_3 &= \frac{-a_1}{2 \cdot 3} = \frac{-a_1}{3!} \\ -n=3: \quad a_5 &= \frac{-a_3}{4 \cdot 5} = \frac{a_1}{5!} \\ -n=5: \quad a_7 &= \frac{-a_5}{6 \cdot 7} = \frac{-a_1}{7!} -\end{align}
-Now, let's put that all together: -\begin{align} -y(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \dots \\ -y &= a_0 \left( 1 - \frac{x^2}{2!} + \frac{x^4}{4!} + \dots \right) + a_1 \left( x - \frac{x^3}{3!} + \frac{x^5}{5!} + \dots \right) -\end{align} -which you might recognize as being the Taylor series expansion of sine and cosine: -\begin{equation} -y(x) = a_0 \cos x + a_1 \sin x -\end{equation} -So, our unknown coefficients end up being our integration constants, which we can use our two constraints to find.
-Find the solution to the ODE -\begin{equation} -(1 + x^2) y^{\prime\prime} - 4x y^{\prime} + 6y = 0 -\end{equation}
-First, rearrange into standard form: -\begin{equation} -y^{\prime\prime} - \frac{4x}{1+x^2} y^{\prime} + \frac{6}{1+x^2} y = 0 -\end{equation} -Then, check whether $x=0$ is an ordinary point: yes, it is.
-Now, let's insert the power series into the ODE: -\begin{align} -y^{\prime\prime} + x^2 y^{\prime\prime} - 4 x y^{\prime} + 6 y &= 0 \\ -\sum_{n=2}^{\infty} a_n (n)(n-1)x^{n-2} + x^2 \sum_{n=0}^{\infty} a_n (n)(n-1)x^{n-2} - 4 x \sum_{n=1} a_n (n) x^{n-1} + 6 \sum_{n=0}^{\infty} a_n x^n &= 0 -\end{align} -First, we'll use the power rule: -\begin{equation} -\sum_{n=2}^{\infty} a_n (n)(n-1)x^{n-2} + \sum_{n=2}^{\infty} a_n (n)(n-1)x^{n} - 4 \sum_{n=1} a_n (n) x^{n} + 6 \sum_{n=0}^{\infty} a_n x^n = 0 -\end{equation} -and then the index shift and dummy index rules on the first term: -\begin{equation} -\sum_{n=2}^{\infty} a_n (n)(n-1)x^{n-2} \rightarrow \sum_{m=0}^{\infty} a_{m+2} (m+2)(m+1) x^m \rightarrow \sum_{n=0}^{\infty} a_{n+2} (n+2)(n+1) x^n -\end{equation}
-Then, put that back into the full equation and combine the sums: -\begin{align} -\sum_{n=0}^{\infty} a_{n+2} (n+2)(n+1) x^n + \sum_{n=0}^{\infty} a_n (n)(n-1)x^{n} - 4 \sum_{n=1} a_n (n) x^{n} + 6 \sum_{n=0}^{\infty} a_n x^n &= 0 \\ -\sum_{n=0}^{\infty} x^n \left[ a_{n+2} (n+2)(n+1) + a_n (n)(n-1) - 4a_n (n) + 6a_n \right] &= 0 \\ -a_{n+2} (n+2)(n+1) + a_n (n^2 -5n + 6) &= 0 \\ -a_{n+2} (n+2)(n+1) + a_n (n-3)(n-2) &= 0 \\ -\end{align}
-Thus, our recursion formula for the coefficients $a_n$ is -\begin{equation} -a_{n+2} = -a_n \frac{(n-3)(n-2)}{(n+1)(n+2)} -\end{equation} -Again, we can see that the even terms will be related and the odd terms will be related: -\begin{align} -n=0: \quad a_2 &= -a_0 \frac{6}{2} = -3 a_0 \\ -n=2: \quad a_4 &= 0 \\ -n=4: \quad a_6 &= -a_4 \frac{2}{30} = 0 \\ -&\ldots -\end{align} -and the odd terms: -\begin{align} -n=1: \quad a_3 &= -a_1 \frac{2}{6} = \frac{-a_1}{3} \\ -n=3: \quad a_5 &= 0 \\ -n=5: \quad a_7 &= -a_5 \frac{6}{42} = 0 \\ -&\ldots -\end{align}
-The solution is then -\begin{align} -y(x) &= a_0 + a_1 x - 3 a_0 x^2 - \frac{a_1}{3} x^3 \\ -y &= a_0 \left(1 - 3x^2 \right) + a_1 \left( x - \frac{x^3}{3} \right) -\end{align} -where we find $a_0$ and $a_1$ using our initial or boundary conditions.
- -This chapter focuses on analytical and numerical methods for solving 2nd-order ordinary differential equations (ODEs), including initial-value problems (IVPs) and boundary-value problems (BVPs).
- -'+result.length+' Result(s) found
'); - for (var item in result) { - var ref = result[item].ref; - if(store[ref].teaser){ - var searchitem = - ''+store[ref].excerpt.split(" ").splice(0,20).join(" ")+'...
'+ - ''+store[ref].excerpt.split(" ").splice(0,20).join(" ")+'...
'+ - '
-
-