Scientific Computing - Quiz 4 - Problem 3-13

Problem:

To be copied

Solution:

The differential equation is so easy so we can first solve it analytically, again, using Laplace transform.

First, let $ y(x) = \psi(x - 1) $, that allow us to pose the initial condition at $ x = -1 $ to $ x = 0 $.

$ \begin{eqnarray*} y'' &=& \psi''(x - 1) \\ &=& (\beta - 100)\psi(x - 1) \\ &=& (\beta - 100)y \\ \mathcal{L}(y'') &=& \mathcal{L}((\beta - 100)y) \\ s\mathcal{L}(y') - y'(0) &=& (\beta - 100)\mathcal{L}(y) \\ s\mathcal{L}(y') - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s(s\mathcal{L}(y) - y(0)) - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s(s\mathcal{L}(y) - 0) - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s^2\mathcal{L}(y) - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s^2\mathcal{L}(y) + (100 - \beta)\mathcal{L}(y) &=& 1 \\ \mathcal{L}(y) &=& \frac{1}{s^2 + (100 - \beta)} \\ \mathcal{L}(y) &=& \frac{1}{\sqrt{(100 - \beta)}}\frac{\sqrt{(100 - \beta)}}{s^2 + (100 - \beta)} \\ y &=& \frac{1}{\sqrt{(100 - \beta)}}\sin(\sqrt{(100 - \beta)}x) \\ \psi(x) &=& \frac{1}{\sqrt{(100 - \beta)}}\sin(\sqrt{(100 - \beta)}(x + 1)) \\ \end{eqnarray*} $

Here is the MATLAB code for solving the ODE:

clear all; close all; clc

main driver (just paste into MATLAB)
tol = 10^(-4);
xspan = [-1 1];
A = 1;
ic = [0 A];
beta = 91;
[t y] = ode45('f_rhs', xspan, ic, [], beta);
plot(t, y(:,1))
y(end, 1)

and the function to be solved:

function rhs=f_rhs(x, ic, dummy, beta)
  y1 = ic(1);
  y2 = ic(2);
  rhs = [y2;
         (beta - 100)*y1];

this is the overall code for finding all eigenvalues

clear all;                                         % clear all previously defined variables
close all;                                         % clear all previously defined figures

tol=10^(-4);                                       % define a tolerance level to be achieved by the shooting algorithm
col=[’r’,’b’,’g’,’c’,’m’,’k’];                     % eigenfunction colors
n0=100;                                            % define the parameter n0
A=1;                                               % define the initial slope at x=-1
x0=[0 A];                                          % initial conditions: x1(-1)=0, x1’(-1)=A
xp=[-1 1];                                         % define the span of the computational domain
beta_start=n0;                                     % beginning value of beta

for modes=1:5                                      % begin mode loop
  beta=beta_start;                                 % initial value of eigenvalue beta
  dbeta=n0/100;                                    % default step size in beta
  for j=1:1000                                     % begin convergence loop for beta
    [t,y]=ode45(’f_rhs’,xp,x0,[],n0,beta);         % solve ODEs
    if abs(y(end,1)-0) < tol                       % check for convergence
      beta                                         % write out eigenvalue
      break                                        % get out of convergence loop
    end
    if (-1)^(modes+1)*y(end,1)>0                   % this IF statement block
      beta=beta-dbeta;                             % checks to see if beta
    else                                           % needs to be higher or lower
      beta=beta+dbeta/2;                           % and uses bisection to
      dbeta=dbeta/2;                               % converge to the solution
    end %
  end % end convergence loop
  beta_start=beta-0.1;                             % after finding eigenvalue, pick new starting value for next mode
  norm=trapz(t,y(:,1).*y(:,1))                     % calculate the normalization
  plot(t,y(:,1)/sqrt(norm),col(modes)); hold on    % plot modes
end                                                % end mode loop

Now we can check with the analytic answers - note the output beta values actually lead to full/half cycle of the sine.

In order for $ \psi(1) = 0 $, we need $ \frac{1}{\sqrt{(100 - \beta)}}\sin(2\sqrt{(100 - \beta)}) = 0 $, of course that correspond to $ \sqrt{100 - \beta} $ is a multiple of $ \frac{\pi}{2} $

Scientific Computing - Quiz 4 - Problem 1, 2

Problem:

Consider the differential equation
$ y''' + y'y +y^2 = \cos(t) $

which is the correct conversion of this equation into first-order form

Consider also the boundary conditions

$ y′′(0)=5,y(0)=3,y(1)y′(1)+y(1)=0 $

What is the correct implementation of these boundary conditions in the first-order system.

Solution:

Let $ y_1 = y $, $ y_2 = y' $, $ y_3 = y'' $, we have

$ y_1' = y_2 $
$ y_2' = y_3 $
$ y_3' = y''' = \cos(t) - y'y - y^2 = \cos(t) - y_1y_2 - y_1^2 $

$ y''(0) = 5 \implies y_3(0) = 5 $
$ y(0) = 3 \implies y_1(0) = 3 $
$ y(1)y′(1)+y(1)=0 \implies y_1(1)y_2(1) +y_1(1) = 0 $

Scientific Computing - Quiz 3 - Question 1-3

Problem:

Solve the simple boundary value problem
$ y′′=−4y $ with $ y(0)=0 $, $ y′(1)=\cos(2) $

(Use your above work and solution to continue with this question)

With boundary value problem considered previously
$ y′′=−4y $ with $ y(0)=0 $, $y′(1)=\cos(2) $

Consider instead

$ y′′=−4y $ with $ y(0) = 0 $, $ y′(0) = A $ and use a bisection to compute the iterative values of $ A_1 = 0.5 $ that converge to the solution using a shooting algorithm (NOTE: use your known solution in the interval $t \in [0,1]$. Here assume that $ A_1 = 0.5 $ (too low) and $ A_2 = 2 $ (too high) to compute (a) $A_3$, (b) $A_4$ and (c) $A_5$.

Solution:

Taking Laplace transform on both side gives

$ \begin{eqnarray*} \mathcal{L}(y'') &=& \mathcal{L}(-4y) \\ s\mathcal{L}(y') - y'(0) &=& -4\mathcal{L}(y) \\ s(s\mathcal{L}(y) - y(0))- y'(0) &=& -4\mathcal{L}(y) \\ s(s\mathcal{L}(y))- y'(0) &=& -4\mathcal{L}(y) \\ s^2\mathcal{L}(y) + 4\mathcal{L}(y) &=& y'(0) \\ \mathcal{L}(y) &=& \frac{y'(0)}{s^2 + 4} \\ \mathcal{L}(y) &=& \frac{y'(0)}{2}\frac{2}{s^2 + 4} \\ y(t) &=& \frac{y'(0)}{2}\sin(2t) \\ \end{eqnarray*} $

Now we see the general form of $ y $, using the boundary condition we can solve for the unknown $ y'(0) $ as follow:

$ \begin{eqnarray*} y'(t) &=& y'(0)\cos(2t) \\ y'(1) &=& y'(0)\cos(2) \\ \cos(2) &=& y'(0)\cos(2) \\ y'(0) &=& 1 \\ y(t) &=& \frac{1}{2}\sin(2t) \\ \end{eqnarray*} $

The rest is just bisection.

Note that $ y'(1) = Acos(2) $

$ A $	$ y'(1) $	goal	$ y'(1) $ is …	$ A $ is …
0.5000	-0.2081	-0.4161	too high	too low
2.0000	-0.8323	-0.4161	too low	too high
1.2500	-0.5202	-0.4161	too low	too high
0.8750	-0.3641	-0.4161	too high	too low
1.0625	-0.4422	-0.4161	too low	too high

Note that the question is confusing me saying $ A = 0.5 $ is too low while $ y'(1) $ when $ A = 0.5 $ is actually too high :(

Also be very careful with the requirement for the number of digits.

Scientific Computing - Quiz 2 - Question 3-6

Problem:

Consider the differential equation
$ \frac{dy}{dt} = y $

with y(0)=1. The exact solution to this problem is given by $ y(t)=e^t $. Knowing the exact solution and using $\Delta t=0.1$ with the Euler scheme, compute the local error for two iterative steps ((a) $\epsilon_1$, (b) $\epsilon_2$) and compute the global error ((c) $E_1$, (d) $E_2$) at each of the three steps.

Solution:

The local error computation is exactly the same as the last problem.

$ \epsilon_1 = e^{0.1} - (1 + 0.1) \approx  0.0052 $
$ \epsilon_2 = e^{0.2} - (e^{0.1} + 0.1e^{0.1}) \approx 0.0057 $

For the 1st step, the global error is actually the same as the local error.

$ E_1 = \epsilon_1 = 0.0052 $

For the 2nd step, we need to actually iterate twice.

$ y_1 = 1 + 0.1 $
$ y_2 = y_1 + 0.1 y_1 = 1.21 $

Therefore $ E_2 = e^{0.2} - 1.21 \approx 0.011 $

Notice an important observation here - $ \epsilon_2 \neq E_2 - E_1 $. The error introduced in the last step also lead to wrong slope estimation, the local error of the second step assumed the correct slope will be there!

Scientific Computing - Quiz 2 - Question 1,2

Problem:

Question 1:

Consider the function $f(t)=\cos(t)$. What is the error bound associated with expanding the function $f(t+\Delta t)$ about $ t $using $\Delta t=0.1$. In particular, compute the error bound when truncated at the second-derivative (keep only the first two significant-digits).

Question 2:

Consider the function $f(t)=\cos(t)$. What is the error bound associated with expanding the function $f(t+\Delta t)$ about $t $ using $\Delta t=0.1$. In particular, compute the error bound when truncated at the third-derivative (keep only the first two significant-digits).

Solution:

Consider the Taylor expansions:

For question 1:

$ \cos(t + \Delta t) = \cos(t) - \Delta t \sin(t) - \frac{\Delta t^2}{2}\cos(c) $
So the error bound is the maximum value of the last term, which is $ \frac{0.1^2}{2} = 0.005 $

For question 2:

$ \cos(t + \Delta t) = \cos(t) - \Delta t \sin(t) - \frac{\Delta t^2}{2}\cos(t) + \frac{\Delta t^3}{6}\sin(c) $
So the error bound is the maximum value of the last term, which is $ \frac{0.1^3}{6} = 0.00017 $

Scientific Computing - Quiz 1 - Question 8

Problem:

Derive the Adams-Moulton scheme and Adams-Bashforth scheme that uses $y_{n+1}, y_{n}, f_{n-1}$ and $y_{n}, f_{n-1}, f_{n-2}$ respectively.

Solution:

To approximate a polynomial with three points, we use the Lagrange's polynomial.

$ \begin{eqnarray*} p(t) &=& f_n\frac{(t - (n-1)\Delta t)(t - (n-2)\Delta t)}{(n \Delta t - (n-1)\Delta t)(n \Delta t - (n-2)\Delta t)} + f_{n-1}\frac{(t - n\Delta t)(t - (n-2)\Delta t)}{((n - 1)\Delta t - n\Delta t)((n - 1)\Delta t - (n-2)\Delta t)} + f_{n-2}\frac{(t - n\Delta t)(t - (n-1)\Delta t)}{((n - 2)\Delta t - n\Delta t)((n - 2)\Delta t - (n-1)\Delta t)} \\ &=& f_n\frac{(t - (n-1)\Delta t)(t - (n-2)\Delta t)}{2 \Delta t^2} + f_{n-1}\frac{(t - n\Delta t)(t - (n-2)\Delta t)}{-\Delta t^2} + f_{n-2}\frac{(t - n\Delta t)(t - (n-1)\Delta t)}{2\Delta t^2} \\ &=& \frac{1}{2\Delta t^2} (f_n(t - (n-1)\Delta t)(t - (n-2)\Delta t) - 2f_{n-1}(t - n\Delta t)(t - (n-2)\Delta t) + f_{n-2}(t - n\Delta t)(t - (n-1)\Delta t)) \\ \end{eqnarray*} $

Note that while the expression is long, it is simply a quadratic polynomial with a lot of awful looking constants. To simplify things, let $ u = t - n\Delta t $

$ \begin{eqnarray*} p(u) &=& \frac{1}{2\Delta t^2} (f_n(u + \Delta t)(u + 2\Delta t) - 2f_{n-1}(u)(u + 2\Delta t) + f_{n-2}(u)(u + \Delta t)) \\ &=& \frac{1}{2\Delta t^2} (f_n(u^2 + 3u\Delta t + 2\Delta t^2) - 2f_{n-1}(u^2 + 2u\Delta t) + f_{n-2}(u^2 + u\Delta t)) \\ &=& \frac{1}{2\Delta t^2} (u^2(f_n - 2f_{n-1} + f_{n-2}) + u\Delta t(3f_n - 4f_{n-1} + f_{n-2}) + 2f_n\Delta t^2) \\ \end{eqnarray*} $

Ugly as is, we can integrate now.

$ \begin{eqnarray*} y_{n + 1} - y_{n} & = & \int\limits_{t_n}^{t_{n+1}}{f(t,y)dt} \\ y_{n + 1} & = & y_{n} + \int\limits_{t_n}^{t_{n+1}}{f(t,y)dt} \\ & \approx & y_{n} + \int\limits_{t_n}^{t_{n+1}}{p(t)dt} \\ & \approx & y_{n} + \int\limits_{0}^{\Delta t}{p(u)du} \\ & \approx & y_{n} + \int\limits_{0}^{\Delta t}{\frac{1}{2\Delta t^2} (u^2(f_n - 2f_{n-1} + f_{n-2}) + u\Delta t(3f_n - 4f_{n-1} + f_{n-2}) + 2f_n\Delta t^2)du} \\ & \approx & y_{n} + \frac{1}{2\Delta t^2} \int\limits_{0}^{\Delta t}{(u^2(f_n - 2f_{n-1} + f_{n-2}) + u\Delta t(3f_n - 4f_{n-1} + f_{n-2}) + 2f_n\Delta t^2)du} \\ & \approx & y_{n} + \frac{1}{2\Delta t^2} (\frac{u^3}{3}(f_n - 2f_{n-1} + f_{n-2}) + \frac{u^2}{2}\Delta t(3f_n - 4f_{n-1} + f_{n-2}) + 2uf_n\Delta t^2) |_{0}^{\Delta t} \\ & \approx & y_{n} + \frac{1}{2\Delta t^2} (\frac{\Delta t^3}{3}(f_n - 2f_{n-1} + f_{n-2}) + \frac{\Delta t^3}{2}(3f_n - 4f_{n-1} + f_{n-2}) + 2\Delta t^3f_n) \\ & \approx & y_{n} + \frac{\Delta t}{2} (\frac{1}{3}(f_n - 2f_{n-1} + f_{n-2}) + \frac{1}{2}(3f_n - 4f_{n-1} + f_{n-2}) + 2f_n) \\ & \approx & y_{n} + \frac{\Delta t}{12} (2(f_n - 2f_{n-1} + f_{n-2}) + 3(3f_n - 4f_{n-1} + f_{n-2}) + 12f_n) \\ & \approx & y_{n} + \frac{\Delta t}{12} (23f_n - 16f_{n-1} + 5f_{n-2}) \\ \end{eqnarray*} $

This is the Adam-Bashforth scheme - I will skip the Adam Moulton scheme, derivation like this is best leave to the computer...

We have the correct answer at https://en.wikipedia.org/wiki/Linear_multistep_method, and it is great to see the answer matches (to be honest to my future self, I did spend great time to fix the derivation to match the solution, the very detailed derivation does help me a lot to fix the bug without redoing everything)

Scientific Computing - Quiz 1 - Question 7

Problem:

Reconsider the fundamental theorem of calculus as expressed in (1.1.16 found in the lecture note packet). But now, integrate from $t-\Delta t$ to $t+\Delta t$. Using (1.1.19), what is the new version of (1.1.20):

Solution:

We start with $ y(t + \Delta t) - y(t - \Delta t) = \int\limits_{-t - \Delta t}^{t + \Delta t}{f(t, y)dt} $.

Next we approximate the integral by constant and put terms around, we get

$ y_{n+1} = y_{n-1} + 2\Delta t f(t, y_n) $.

Scientific Computing - Quiz 1 - Question 6

Problem:

Walk through one time-step of 4th-order Runge-Kutta with simple differential equation system $y'=y$, $y(0)=1$ and $\Delta t=2$.
Identify the correct values for f1; f2; f3; f4 and y(2).

Solution:

It is always good to know the solution to the differential equation first, in this case it is really easy.

$ \begin{eqnarray*} \frac{dy}{dt} &=& y \\ \frac{1}{y}\frac{dy}{dt} &=& 1 \\ \int{\frac{1}{y}\frac{dy}{dt}dt} &=& \int{1dt} \\ \log y &=& t + C_1 \\ y &=& e^{t + C_1} \\ &=& C_2e^t \\ \end{eqnarray*} $

Now the initial condition decides $ C_2 = 1 $, making the solution simply $ y = e^t $.

The analytic solution is not required, it is just that it is a perfect guide for us to evaluate whether our solution make sense. Now, looking up the Runge-Kutta equations, we have

$ f_1 = f(t_n, y_n) $
$ f_2 = f(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} f_1) $
$ f_3 = f(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} f_2) $
$ f_4 = f(t_n + \Delta t, y_n + \Delta t f_3) $

Next, we substitute $ f(t, y) = y $ in the expressions above, then we get:

$ f_1 = y_n $
$ f_2 = y_n + \frac{\Delta t}{2} f_1 $
$ f_3 = y_n + \frac{\Delta t}{2} f_2 $
$ f_4 = y_n + \Delta t f_3 $

Excel will just happy to crank those answers for us now.

The final answer is

$ f_1 = 1 $
$ f_2 = 2 $
$ f_3 = 3 $
$ f_4 = 7 $
$ y(2) = 7 $.

With excel, I fiddled around the parameter. In particular, I tried to tune $ \Delta t $ to 0.1 and iterate 20 times, the result has $ 1.1 \times 10^{-5} $ error only. It is in fact quite accurate.

As an aside to the analytic solution - the integral of $ \int{\frac{1}{y}dy} = \log{y} $ only make sense when $ y > 0 $. To make sense of cases when $ y < 0 $, we can take the sign out first before the integration and the expression become $ \log|y| = sign(y)t + C_1 \implies e^{|y|} = e^{|sign(y)|}e^{t + C_1} \implies e^{y} = e^{t + C_1} $ and the rest follows as is.

Scientific Computing - Quiz 1 - Question 5

Problem:

Using Eqs.~(1.1.10)-(1.1.12) found in the lecture notes found under the course resource link, derive the three important schemes: Heun's method (A = 1/2), modified Euler-Cauchy (A = 0) and Euler (A = 1).

Solution:

Here is the link to the full lecture notes:

The derivation based on 1.1.12 is pretty easy, for Heun's method:

$ A + B = 1 \implies B = \frac{1}{2} $.
$ PB = \frac{1}{2} \implies P = 1 $.
$ BQ = \frac{1}{2} \implies Q = 1 $.

For Euler-Cauchy method:

$ A + B = 1 \implies B = 1 $.
$ PB = \frac{1}{2} \implies P = \frac{1}{2} $.
$ BQ = \frac{1}{2} \implies Q = \frac{1}{2} $.

For Euler method:

$ A + B = 1 \implies B = 0 $.
P and Q are meaningless in the context of the iteration, we do not want any correction at all in this case.

Scientific Computing - Quiz 1 - Question 4

Problem:

Using a Taylor expansion for $ y(t+\Delta t) = y_{n+1} $ and $ y(t-\Delta t)=y_{n-1} $, derive an Euler-type formula for the future state of the system $ y' = f(y,t) $.

Solution:

The Taylor expansions are:

$ y_{n+1} = y(t + \Delta t) = y(t) + \Delta t \frac{dy}{dt}(t) + \frac{\Delta t^2}{2} \frac{d^2y}{dt^2}(c_1) = y(t) + \Delta t f(y, t) + \frac{\Delta t^2}{2} \frac{d^2y}{dt^2}(c_1) $
$ y_{n-1} = y(t - \Delta t) = y(t) - \Delta t \frac{dy}{dt}(t) + \frac{\Delta t^2}{2} \frac{d^2y}{dt^2}(c_2) = y(t) - \Delta t f(y, t) + \frac{\Delta t^2}{2} \frac{d^2y}{dt^2}(c_2) $

Now subtract the equations to get our answer
$ \begin{eqnarray*} y_{n+1} - y_{n-1} &=& (y(t) + \Delta t f(y, t) + \frac{\Delta t^2}{2} \frac{d^2y}{dt^2}(c_1)) - (y(t) - \Delta t f(y, t) + \frac{\Delta t^2}{2} \frac{d^2y}{dt^2}(c_2)) \\ &=& 2\Delta t f(y, t) + \frac{\Delta t^2}{2}(\frac{d^2y}{dt^2}(c_1) -\frac{d^2y}{dt^2}(c_2)) \\ y_{n+1} &=& y_{n-1} + 2\Delta t f(y, t) + \frac{\Delta t^2}{2}(\frac{d^2y}{dt^2}(c_1) -\frac{d^2y}{dt^2}(c_2)) \\ & \approx & y_{n-1} + 2\Delta t f(y, t) \end{eqnarray*} $

Scientific Computing - Quiz 1 - Question 3

Problem:

L'Hopital's rule is a method for determining the value of indeterminate forms. Determine the value of the following limits
$ \lim\limits_{x \to 0} \frac{e^x - 1 - x}{x^2} $

Solution:

$ \begin{eqnarray*} & & \lim\limits_{x \to 0} \frac{e^x - 1 - x}{x^2} \\ &=& \lim\limits_{x \to 0} \frac{e^x - 1}{2x} \\ &=& \lim\limits_{x \to 0} \frac{e^x}{2} \\ &=& \frac{1}{2} \\ \end{eqnarray*} $

Again, this is easy to validate numerically:

x	f(x)	$ e^2 $
0.1	0.517091808	0.000292
0.01	0.501670842	2.79E-06
0.001	0.500166708	2.78E-08
0.0001	0.500016671	2.78E-10
0.00001	0.500000696	4.85E-13

Scientific Computing - Quiz 1 - Question 2

Problem:

L'Hopital's rule is a method for determining the value of indeterminate forms. Determine the value of the following limits:

$ \lim\limits_{x \to 0} \frac{2\sin x - \sin 2x}{x-\sin x} $

Solution:

$ \begin{eqnarray*} & & \lim\limits_{x \to 0} \frac{2\sin x - \sin 2x}{x - \sin x} \\ &=& \lim\limits_{x \to 0} \frac{2\cos x - 2 \cos 2x}{1 -\cos x} \\ &=& \lim\limits_{x \to 0} \frac{-2\sin x + 4 \sin 2x}{\sin x} \\ &=& \lim\limits_{x \to 0} \frac{-2\sin x + 8 \sin x \cos x}{\sin x} \\ &=& \lim\limits_{x \to 0} -2 + 8 \cos x \\ &=& 6 \end{eqnarray*} $

One can easily verify this numerically:

x	f(x)	$e^2$
0.1	5.988008283	0.000144
0.01	5.999880001	1.44E-08
0.001	5.9999988	1.44E-12
0.0001	6.000000163	2.64E-14
0.00001	5.999979671	4.13E-10

Scientific Computing - Quiz 1 - Question 1

Problem:

L'Hopital's rule is a method for determining the value of indeterminate forms. Determine the value of the following limits

$ \lim\limits_{x \to 0}{\frac{\sin x}{x}} = 1 $

Solution:

It is really funny that the question gives out the answer - the problem already give us the limit!

Alright, let's derive it anyway, if we assume the derivative of sine is cosine, then that is easy. We have:

$ \begin{eqnarray*} & & \lim\limits_{x \to 0}{\frac{\sin x}{x}} \\ &=& \lim\limits_{x \to 0}{\frac{\cos x}{1}} \\ &=& 1 \end{eqnarray*} $

Note that we carefully said that we assume the derivative of sine is cosine above. It is a fact! Why do we bother to carefully assume it? The problem with the above is that it is a logically wrong proof! We actually need this limit to prove the derivative of sine and cosine, it is a circular reasoning.

The correct way of proving the identity should be found in any high school text book. This is based on the inequality $ \theta \ge \sin\theta \ge \theta\cos\theta $. The inequality can be argued geometrically, and then we divide through $ \theta $ and apply squeezing principle to get to the limit we wanted.