Andrew's Exercise Solutions: November 2015

Sunday, November 29, 2015

Integral (1)

Let's do some integrals! The good thing about integral is

(a) It requires some tricks to get it done, an interesting intellectual challenge.
(b) It is easy to check if I have got the right answer, even without the model answer, just differentiate it.

Problem:

$\int{\frac{dx}{\sqrt{x} + \sqrt{x + 1}}}$

Solution:

The trick is rationalizing:

$\begin{eqnarray*} & & \int{\frac{dx}{\sqrt{x} + \sqrt{x + 1}}} \\ &=& \int{\frac{dx}{\sqrt{x} + \sqrt{x + 1}}\frac{\sqrt{x + 1} - \sqrt{x}}{\sqrt{x + 1} - \sqrt{x}}} \\ &=& \int{(\sqrt{x + 1} - \sqrt{x})dx} \\ &=& \frac{2}{3}(x+1)^{\frac{3}{2}} - \frac{2}{3}x^{\frac{3}{2}} + C \end{eqnarray*}$

Saturday, November 28, 2015

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}$

Friday, November 27, 2015

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.

Wednesday, November 25, 2015

Algorithms: Design and Analysis, Part 1 - Index

Algorithms: Design and Analysis, Part 1 - Exam Question 20

Problem:

Recall the Master Method and its three parameters

$a,b,d$ . Which of the following is the best interpretation of

$b^d$ , in the context of divide-and-conquer algorithms?

Solution:

I revisited the lecture video for this one, and now I am explaining in my own words.

The current amount of work is

$O(n^d)$ .
The problem size become

$\frac{n}{b}$ , and then the amount of work (for the merge step, ignoring recursive calls is

$O((\frac{n}{b})^d)$ .

Therefore the amount of work is reduced by a factor of

$O(b^d)$ , that's why this option is correct.

(Option 1) The rate at which the work-per-subproblem is shrinking (per level of recursion).

Algorithms: Design and Analysis, Part 1 - Exam Question 19

Problem:

Running time of Strassen's matrix multiplication algorithm: Suppose that the running time of an algorithm is governed by the recurrence

$T(n)=7*T(n/2)+n$ 2

$. What's the overall asymptotic running time (i.e., the value of$ T(n) $)?

Solution:

Applying master theorem, we have the variables

$a = 7, b = 2, d = 2$ , and we have the

$7 = a > b^d = 4$ case, therefore the answer is

$\theta(n^{\log_2{7}})$

Algorithms: Design and Analysis, Part 1 - Exam Question 18

Problem:

Suppose you are given

$k$ sorted arrays, each with

$n$ elements, and you want to combine them into a single array of kn elements. Consider the following approach. Divide the

$k$ arrays into

$k/2$ pairs of arrays, and use the Merge subroutine taught in the MergeSort lectures to combine each pair. Now you are left with

$k/2$ sorted arrays, each with

$2n$ elements. Repeat this approach until you have a single sorted array with kn elements. What is the running time of this procedure, as a function of

$k$ and

$n$ ?

Solution:

The first iteration deal with

$\frac{k}{2}$ array of length

$n$ , so it takes

$\theta(nk)$ time.
The second iteration deal with

$\frac{k}{4}$ array of length

$2n$ , so it also takes

$\theta(nk)$ time.

Therefore the overall execution time is really just

$\theta(nk\log k)$ time.

This is simply and I should have got this correct, I made a careless mistake here, sigh.

Algorithms: Design and Analysis, Part 1 - Exam Question 17

Problem:

Which of the following statements about Dijkstra's shortest-path algorithm are true for input graphs that might have some negative edge lengths?

Solution:

(Option 1) It is guaranteed to terminate.
The algorithm, regardless of correctness, also make progress towards completion by removing a vertex from the unknown set, so it must terminate.

(Option 2) It may or may not correctly compute shortest-path distances (from a given source vertex to all other vertices), depending on the graph.
This is true because the graph may not have any negative edge at all.

(Option 3) It may or may not terminate (depending on the graph).
This is not true given option 1

(Option 4) It is guaranteed to correctly compute shortest-path distances (from a given source vertex to all other vertices).
This is not true, or else we don't really need the restriction.

Algorithms: Design and Analysis, Part 1 - Exam Question 16

Problem:

Which of the following patterns in a computer program suggests that a hash table could provide a significant speed-up (check all that apply)?

Solution:

Hash table are great for lookup but do not good for finding minimum of maximum.

Algorithms: Design and Analysis, Part 1 - Exam Question 15

Problem:

Which of the following patterns in a computer program suggests that a heap data structure could provide a significant speed-up (check all that apply)?

Solution:

Heap is good for finding minimum and maximum but do no good to lookups

Algorithms: Design and Analysis, Part 1 - Exam Question 14

Problem:

Suppose you implement the operations Insert and Extract-Min using a sorted array (from biggest to smallest). What is the worst-case running time of Insert and Extract-Min, respectively? (Assume that you have a large enough array to accommodate the Insertions that you face.)

Solution:

This is again a done assignment problem, find the solution below

Algorithms: Design and Analysis, Part 1 - Problem Set 5 - Question 3

Algorithms: Design and Analysis, Part 1 - Exam Question 13

Problem:

Suppose that a randomized algorithm succeeds (e.g., correctly computes the minimum cut of a graph) with probability

$p$ (with

$0 < p < 1$ ). Let

$\epsilon$ be a small positive number (less than 1). How many independent times do you need to run the algorithm to ensure that, with probability at least

$1 - \epsilon$ , at least one trial succeeds?

Solution:

Suppose we run the algorithm for

$n$ times and therefore the probability for all of these independent trials to fail is

$(1 - p)^n$ .

The probability to success is therefore

$1 - (1 - p)^n > 1 -\epsilon$ .

$\begin{eqnarray*} 1 - (1 - p)^n &>& 1 -\epsilon \\ - (1 - p)^n &>& -\epsilon \\ (1 - p)^n &<& \epsilon \\ \log((1 - p)^n) &<& \log \epsilon \\ n\log(1 - p) &<& \log \epsilon \\ n &>& \frac{\log \epsilon}{\log(1 - p)} \end{eqnarray*}$

Note the change of inequality sign in the last line, this is because

$\log(1 - p) < 0$

Algorithms: Design and Analysis, Part 1 - Exam Question 12

Problem:

When does a directed graph have a unique topological ordering?

Solution:

(Option 1) Whenever it is directed acyclic
This is wrong 1->2, 3 has two topological order

(Option 2) Whenever it has a unique cycle
This is wrong, any graph with a cycle do not have a topological order.

(Option 3) Whenever it is a complete directed graph
A complete directed graph of size > 1 will necessarily have cycle, so no topological order at all.

(Option 4) None of the other options
This is the only possible solution now, and that's correct!

Algorithms: Design and Analysis, Part 1 - Exam Question 11

Problem:

Which of the following statements hold? (As usual n and m denote the number of vertices and edges, respectively, of a graph.) [Check all that apply.]

Solution:

Breadth-first search can be used to compute the connected components of an undirected graph in

$O(m+n)$ time.
Depth-first search can be used to compute the strongly connected components of a directed graph in

$O(m+n)$ time.
Breadth-first search can be used to compute shortest paths in

$O(m+n)$ time (when every edge has unit length).
Depth-first search can be used to compute a topological ordering of a directed acyclic graph in

$O(m+n)$ time.

All of these options are correct and all of them are just standard results.

Algorithms: Design and Analysis, Part 1 - Exam Question 10

Problem:

Let

$0< \alpha <.5$ be some constant. Consider running the Partition subroutine on an array with no duplicate elements and with the pivot element chosen uniformly at random (as in QuickSort and RSelect). What is the probability that, after partitioning, both subarrays (elements to the left of the pivot, and elements to the right of the pivot) have size at least

$\alpha$ times that of the original array?

Solution:

It is easiest to look at problem of this sort in term of concrete numbers, it is very easy that way, let say

$\alpha = 0.2$ and

$n = 10$ . Then the requirement is that we have at least 2 elements in both subarrays. Further suppose the numbers are simply

$1,2,\cdots, 10$ , then as long as we are not picking 1, 2, 9, 10, we are good to go.

Generalizing, it is easy to see that if we are not picking the

$n\alpha$ largest numbers and

$n\alpha$ smallest numbers, then we satisfy the condition, the probability of that is simply

$1 - 2\alpha$ .

Algorithms: Design and Analysis, Part 1 - Exam Question 9

Problem:

Let f and g be two increasing functions, defined on the natural numbers, with

$f(1)$ ,

$g(1) \ge 1$ . Assume that

$f(n)=O(g(n))$ . Is

$2^{f(n)}=O(2^{g(n)})$ ? (Multiple answers may be correct, check all that apply.)

Solution:

This is simply an assignment problem I have done before, it can be found in the link below.

Algorithms: Design and Analysis, Part 1 - Problem Set 1 - Question 3

The key idea is just going back to the definition and argue the necessary condition for

$2^{f(n)}=O(2^{g(n)})$ is

$f(n) \le g(n)$ .

Algorithms: Design and Analysis, Part 1 - Exam Question 8

Problem:

What is the asymptotic running time of Randomized QuickSort on arrays of length

$n$ , in expectation (over the choice of random pivots) and in the worst case, respectively?

Solution:

It is again just memorizing standard results:

The expected result is

$\theta(n \log n)$ , and the worst case result is

$\theta(n^2)$ .

Algorithms: Design and Analysis, Part 1 - Exam Question 7

Problem:

On adding one extra edge to a directed graph G, the number of strongly connected components...?

Solution:

Anything could happen!

Suppose the graph is a direct chain, adding an edge to make it a cycle can make it reduce a lot of strongly connected components.

1 -> 2 -> 3 -> 4 (4 strongly connected components)

Adding 4 -> 1 will result in just 1 strongly connected components.

On the other hand adding one more edge, say 1 -> 3, does not decrease it any more.

So the answer is this choice:

...might or might not remain the same (depending on the graph).

Algorithms: Design and Analysis, Part 1 - Exam Question 6

Problem:

What is the asymptotic running time of the Insert and Extract-Min operations, respectively, for a heap with

$n$ objects?

Solution:

This is just memorizing the standard results.

Insert takes

$O(\log n)$ time.
Extract-Min takes

$O(\log n)$ time.

Note that the problem is a little sloppy here, it requires an answer in

$\Theta$ , however, in the best case both operation takes just

$O(1)$ time, and the problem do not explicitly ask for the worst case time.

Algorithms: Design and Analysis, Part 1 - Exam Question 5

Problem:

What is the running time of depth - first search, as a function of n and m, if the input graph

$G = (V,E)$ is represented by an adjacency matrix (i.e., NOT an adjacency list), where as usual

$n=|V|$ and

$m=|E|$ ?

Solution:

Each node need to be gone through. And for each of them it need to look at all other nodes for a possible edge, so the time is

$\theta(n^2)$ .

Algorithms: Design and Analysis, Part 1 - Exam Question 4

Problem:

Consider a directed graph

$G = (V,E)$ with non-negative edge lengths and two distinct vertices

$s$ and

$t$ of

$V$ . Let

$P$ denote a shortest path from

$s$ to

$t$ in

$G$ . If we add 10 to the length of every edge in the graph, then:

Solution:

Let consider this graph:

s -1-> a
a -1-> t
s -3-> t

The current shortest path is s -1-> a -1-> t has cost 2.

Suppose we add 10 to each edges, the path is no longer optimal, the optimal path will become

s -13-> t

The general transformation of path cost is path cost + number of edges times 10, therefore, path with less number of edge will be proportionally shorten than the one with larger number of edges, the construction of the example above is based on this.

Therefore we can consider the choices now:

(Option 1) $P$ definitely remains a shortest $s - t$ path.
This is obviously wrong given the example above.

(Option 2) $P$ might or might not remain a shortest $s − t$ path (depending on the graph).
This is true, the example above shown this could be not the shortest path. It could be the shortest path, for example, in case there is only one path.

(Option 3) $P$ definitely does not remain a shortest $s - t$ path.
This is obviously wrong given the above.

(Option 4) If $P$ has only one edge, then $P$ definitely remains a shortest $s − t$ path.
This is true, but it is non-obvious. This can be proved by looking at the cost transformation.

Suppose the contrary that there exists another path

$Q$ after the transformation is shorter than

$P$ , then

New Cost of Q < New Cost of P
Old Cost of Q + num edges in Q times 10 < Old Cost of P + 10

Now num edges in Q must be at least 1, so we have

Old Cost of Q + 10 < Old Cost of P + 10
Old Cost of Q < Old Cost of P

But that is impossible!

$P$ was the shortest path!

Therefore we proved option 4 is true.

Algorithms: Design and Analysis, Part 1 - Exam Question 3

Problem:

What is the asymptotic worst-case running time of MergeSort, as a function of the input array length

$n$ ?

Solution:

Actually all recursion levels will take

$\theta(n)$ time to create the merged arrays, so the time is obviously

$\theta(n \log n)$ time.

Algorithms: Design and Analysis, Part 1 - Exam Question 2

Problem:

Here is an array of ten integers: 5 3 8 9 1 7 0 2 6 4

Suppose we run MergeSort on this array. What is the number in the 7th position of the partially sorted array after the outermost two recursive calls have completed (i.e., just before the very last Merge step)? (When we say "7th" position, we're counting positions starting at 1; for example, the input array has a "0" in its 7th position.)

Solution:

Just before the last merge, the left half and the right half must be sorted, so it should be

1 3 5 8 9 0 2 4 6 7

Therefore the solution is 2.

Algorithms: Design and Analysis, Part 1 - Exam Question 1

Problem:

Recall the Partition subroutine that we used in both QuickSort and RSelect. Suppose that the following array has just been partitioned around some pivot element: 3, 1, 2, 4, 5, 8, 7, 6, 9

Solution:

Any number that can be a pivot must have all left hand side smaller than it and right hand side larger than it, looking at the array, that will be these elements:

3, 1, 2, 4, 5, 8, 7, 6, 9

Therefore the answer is 4, 5, 9. I missed 9 during my own exam, just a careless mistake.

Monday, November 23, 2015

Focus and directix for a parabola

A parabola is defined as the locus of point where its distance to a line (called directrix) and a point (called focus) are equal.

As an canonical example, we define the directrix to be

$x = -a$ , and the focus be

$(a, 0)$ , then the parabola will be defined as point such that

$(x + a)^2 = (x - a)^2 + y^2$

Simplifying we will get

$y^2 = 4ax$ .

The interesting thing is to find directix and focus from an arbitrary (let's say, still axis parallel) parabolas. It is best illustrated with an example

Suppose the parabola is

$-x^2 - 2y + 2x = 12$

Now we have

$-x^2$ , so the parabola open downwards, it is best to deal with the positive coefficients, so we do:

$-12 - 2y = x^2 - 2x$

Next we complete the square, this is trivial

$-11 - 2y = (x - 1)^2$

Last we just fit it into the form we want:

$4\frac{1}{2}(y + \frac{11}{2}) = -(x - 1)^2$

So we know it is just a shift of the standard form

$4ay = x^2$ , the vertex is

$(1, \frac{-11}{2})$ , the focus is

$(1, -6)$ , the directrix is

$y = -5$ . The focus and directix are found by moving

$a = \frac{1}{2}$ units from the vertex.

Here is a plot showing the elements

Last but not least, it is always nice to check answer:

$(y + 5)^2 = (x - 1)^2 + (y + 6)^2$

$y^2 + 10y + 25 = x^2 - 2x + 1 + y^2 + 12y + 36$

$-x^2 - 2y + 2x = 12$

So yes, we have got the right answer!

Saturday, November 21, 2015

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

Friday, November 20, 2015

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.

Algorithms: Design and Analysis, Part 1 - Programming Question 6 - Part 2

Problem:

Keep track of the median when data arrives 1 by 1, compute the sum of all medians modulo 10000 for each element.

Solution:

I loved this idea in the lecture, we maintain two heaps. The low heap contains the smaller half of the integers and we enable extract max there. The high heap contains the bigger half of the integers and we enable extract min there. We can use constant number of heap operation to maintain the invariant and therefore the per data item cost is

$O(\log n)$ .

	#include <iostream>
	#include <fstream>
	#include <algorithm>

	using namespace std;

	const char* filename = "median.txt";
	const int NUM_ELEMENTS = 10000;

	class heap
	{
	public:
	heap();
	int get_min();
	void delete_min();
	void insert(int x);
	void fixup(int index);
	void fixdown(int index);
	private:
	int data[NUM_ELEMENTS];
	int size;
	};

	heap::heap()
	{
	size = 0;
	}

	void heap::delete_min()
	{
	if (size == 0)
	{
	throw 0;
	}
	else
	{
	data[1] = data[size];
	size--;
	fixdown(1);
	}
	}

	int heap::get_min()
	{
	if (size == 0)
	{
	throw 0;
	}
	else
	{
	return data[1];
	}
	}

	void heap::insert(int x)
	{
	size++;
	data[size] = x;
	fixup(size);
	}

	void heap::fixup(int index)
	{
	if (index == 1)
	{
	return;
	}
	else
	{
	int parent = index / 2;
	if (data[parent] > data[index])
	{
	swap(data[parent], data[index]);
	fixup(parent);
	}
	}
	}

	void heap::fixdown(int index)
	{
	int left_child = index * 2;
	int right_child = left_child + 1;
	if (right_child <= size)
	{
	if (data[left_child] < data[right_child])
	{
	if (data[index] < data[left_child])
	{
	// we are good, index < left < right
	}
	else
	{
	swap(data[index], data[left_child]);
	fixdown(left_child);
	}
	}
	else
	{
	if (data[index] < data[right_child])
	{
	// we are good, index < right < left
	}
	else
	{
	swap(data[index], data[right_child]);
	fixdown(right_child);
	}
	}
	}
	else if (left_child <= size)
	{
	if (data[index] < data[left_child])
	{
	// we are good, index < left < right
	}
	else
	{
	swap(data[index], data[left_child]);
	fixdown(left_child);
	}
	}
	}

	int main()
	{
	int median = 0;
	heap low;
	heap high;
	fstream fin;
	fin.open(filename);

	int median_sum = 0;

	for (int i = 1; i <= NUM_ELEMENTS; i++)
	{
	int data;
	fin >> data;
	if (i == 1)
	{
	median = data;
	}
	else if (i == 2)
	{
	if (data < median)
	{
	low.insert(-data);
	high.insert(median);
	}
	else
	{
	low.insert(-median);
	high.insert(data);
	}
	}
	else if (i % 2 == 1)
	{
	int low_max = -low.get_min();
	int high_min = high.get_min();

	if (data < low_max)
	{
	low.delete_min();
	low.insert(-data);
	median = low_max;
	}
	else if (data > high_min)
	{
	high.delete_min();
	high.insert(data);
	median = high_min;
	}
	else
	{
	median = data;
	}
	}
	else if (i % 2 == 0)
	{
	int low_max = -low.get_min();
	int high_min = high.get_min();

	if (data < low_max)
	{
	low.insert(-data);
	high.insert(median);
	}
	else if (data > high_min)
	{
	high.insert(data);
	low.insert(-median);
	}
	else
	{
	if (data > median)
	{
	high.insert(data);
	low.insert(-median);
	}
	else
	{
	high.insert(median);
	low.insert(-data);
	}
	}
	}

	if (i % 2 == 0)
	{
	median = -low.get_min();
	}

	median_sum = (median_sum + median) % 10000;
	}
	cout << median_sum << endl;
	fin.close();
	return 0;
	}

view raw Main.cpp hosted with ❤ by GitHub

Algorithms: Design and Analysis, Part 1 - Programming Question 6 - Part 1

Problem:

Find the number of integers in [-10,000, 10,000] such that it is representable as a sum of two distinct integers in a list of a million integers.

Solution:

The original solution is to use hash table to look up, but that is too slow for this problem, people reported that will take hours. The distribution of the number has a very wide range and therefore sorting helped a lot, in particular, we will keep two pointers, left and right

If data[left] + data[right] > 10,000, there is no way there are any solution starting with data[right], so we decrement right by 1.

Similarly, if data[left] + data[right] < -10,000, we can increment left by 1.

Otherwise we found a match, but then what - we cannot discard left/right.

Now just move on a try if data[start] can match anything on the right by scanning right inward without moving the actual right pointer yet. Similarly we do the same for data[end].

With the above exercise, we should be exhausted everything that can start with start or end, so we are done with them, we decrement right by 1 and increment left by 1.

Here is the code of that!

	#include <algorithm>
	#include <iostream>
	#include <fstream>

	using namespace std;

	typedef long long int64;

	const char* filename = "algo1-programming_prob-2sum.txt";
	const int NUM_ELEMENTS = 1000000;
	const int MIN = -10000;
	const int MAX = 10000;

	int64 data[NUM_ELEMENTS];

	int main()
	{
	ifstream fin;
	fin.open(filename);

	for (int i = 0; i < NUM_ELEMENTS; i++)
	{
	fin >> data[i];
	}

	sort(data, data + NUM_ELEMENTS);

	int start = 0;
	int end = NUM_ELEMENTS - 1;
	bool found[MAX - MIN + 1];
	for (int i = MIN; i <= MAX; i++)
	{
	found[i - MIN] = false;
	}
	while (start < end)
	{
	int64 probe_sum = data[start] + data[end];
	if (probe_sum < MIN)
	{
	// the value is too small, there is just no hope for success, let go the small side
	start++;
	}
	else if (probe_sum > MAX)
	{
	// the value is too large, there is just no hope for success, let go the small side
	end--;
	}
	else
	{
	if (data[start] != data[end])
	{
	found[probe_sum - MIN] = true;
	}
	int current_start = start;
	int current_end = end;
	while (true)
	{
	// let see if there are any more solution starting with the same end
	start++;
	int64 probe_sum = data[start] + data[end];
	if (probe_sum < MIN)
	{
	// This is impossible
	break;
	}
	else if (probe_sum > MAX)
	{
	break;
	}
	else
	{
	if (data[start] != data[end])
	{
	found[probe_sum - MIN] = true;
	}
	}
	}
	start = current_start;

	while (true)
	{
	// let see if there are any more solution starting with the same start
	end--;
	int64 probe_sum = data[start] + data[end];
	if (probe_sum < MIN)
	{
	break;
	}
	else if (probe_sum > MAX)
	{
	// This is impossible
	break;
	}
	else
	{
	if (data[start] != data[end])
	{
	found[probe_sum - MIN] = true;
	}
	}
	}
	end = current_end;
	// We have exhausted all solution with start and end, so skip them both
	start++;
	end--;
	}
	}

	int count = 0;
	for (int i = MIN; i <= MAX; i++)
	{
	if (found[i - MIN])
	{
	count++;
	}
	}
	cout << count << endl;

	fin.close();
	return 0;
	}

view raw Main.cpp hosted with ❤ by GitHub