assignment_6
MSE3114 Computational Methods for Physicists and Materials Engineers
Assignment 6
Q1. (Understanding algorithm by graphics) Show the procedures to find fixed point or root by graphics. Draw the lines by your hands directly on the following figures, rather than by computer programming.
Given the initial guess $x_0$, find $x_1$, $x_2$, $x_3$ and $x_4$ following the iteration algorithm of finding fixed point: $x_{n+1} = f(x_n)$. According to your drawing, will this iteration converge?
Given the initial guess $x_0$, find $x_1$, $x_2$, $x_3$ and $x_4$ following the iteration algorithm of finding fixed point: $x_{n+1} = f(x_n)$. According to your drawing, will this iteration converge?
Given the initial guess $x_0$, find $x_1$, $x_2$ and $x_3$ following the Newton’s algorithm of root finding. Give a brief explanation to how you draw the lines and read $\{x_n\}$ from the figure. According to your drawing, will this iteration converge?
Given the initial guess $x_0$, find $x_1$, $x_2$ and $x_3$ following the chord method (quasi-Newton) of root finding. Give a brief explanation to how you draw the lines and read $\{x_n\}$ from the figure. According to your drawing, will this iteration converge?
(Bonus)
Given the initial guess $x_0$, find $x_1$, $x_2$ and $x_3$ following the secant method (quasi-Newton) of root finding. Give a brief explanation to how you draw the lines and read $\{x_n\}$ from the figure. According to your drawing, will this iteration converge?
(10 + 10 + 10 + 10 + (bonus 5) marks)
Q2. (Root finding) Solve the following equations by python codes.
Consider the nonlinear equation:
$$
x^4 + 7 x^3 + 2 x^2 + x + 1 = 0.
$$
Firstly, plot the functions $f_1(x) = x^4 + 7 x^3 + 2 x^2 + x + 1$ and $f_2(x) = 0$ by matplotlib. The suggested range for $x$ is from $-7$ to $2$. From the plot, find the approximate roots. Secondly, try different initial guesses and find all the roots to the equation by scipy.optimize.root(). You can check the solutions with the approximate roots observed from the plot.
Consider the nonlinear equation:
$$
\sin x = \frac{1}{x}.
$$
Firstly, plot the functions $f_1(x) = \sin x$ and $f_2(x) = 1/x$ by matplotlib. The suggested range for $x$ is from $0.1$ to $20$. From the plot, find the approximate roots. Secondly, try different initial guesses and find all the roots to equation in the range $x \in [0.1, 20]$ by scipy.optimize.root(). You can check the solutions with the approximate roots observed from the plot.
Consider the system of nonlinear equations:
\begin{equation}
\left\{\begin{array}{l}
xy = 10 \\
x^2 + y^2 = 36
\end{array}\right..
\end{equation}
Firstly, plot the two functions by matplotlib. The suggested range is $x \in [-7, 7]$, $y \in [-7, 7]$. From the plot, find the approximate roots. Secondly, try different initial guesses and find all the roots to the equation by scipy.optimize.root(). You can check the solutions with the approximate roots observed from the plot.
(10 + 10 + 10 marks)
Q3. (Fixed point) Johannes Kepler calculated the orbit of a planet around the sun by solving a nonlinear equation, i.e., Kepler’s equation:
$$
f(E) \equiv M + e \sin E = E,
$$
where $M$ is the mean anomaly, $E$ is the eccentric anomaly and $e$ is the eccentricity.
Solving this equaiton for $E$ is a fixed point problem.
In Assignment 5, we have proved that when $|e| < 1$, $f$ is a contraction and a fixed point exists.
Assume that the orbit of a planet around the sun is elliptic.
For this elliptic orbit, the semi-major axis is $a$ and the eccentricity is $e$.
The period, i.e., the time for the planet to complete one cycle, is $P$.
Then, at a moment $t$, $M(t) = 2\pi t/P$.
$E(t)$ can be obtained by finding the fixed point for the Kepler's equation.
Finally, at the moment $t$, the distance between the planet and the sun is
$$
r = a (1 - e\cos E).
$$
Following the above procedure, write a python code to calculate the distance between the earth and the sun as a function of time. Use matplatlib to plot this function. The parameters are: $e = 0.0167$, $a = 1$ AU (Astronomical Unit), and $P = 365.2564$ days.
(15 marks)
Q4. (Root finding) The Schrödinger equation for a single particle in one dimensional space is
$$
-\frac{\hbar^2}{2m} \frac{d^2 \psi(x)}{d x^2} + V_0 \psi(x) = E \psi(x),
$$
where $\hbar$ is the reduced Planck constant, $m$ is the mass, $\psi(x)$ is the wave function, $V_0$ is the potential energy and $E$ is the energy.
The solution to this equation depends on the boundary conditions.
Now, we consider a particle in a finite potential well as illustrated in the following figure.
Finally, we find that the energy satisfying either of the following nonlinear equations is the energy allowed by the system.
(i) The energy for symmetric wave function satisfies
$$
\sqrt{E} \tan\left(\frac{\pi}{2}\sqrt{E}\right) = \sqrt{V_0 - E}. \quad(1)
$$
(ii) The energy for asymmetric wave function satisfies
$$
-\sqrt{E} \cot\left(\frac{\pi}{2}\sqrt{E}\right) = \sqrt{V_0 - E}. \quad(2)
$$
Since the solutions to either equation are discretized set of values, the energy allowed by the system is quantized.
The questions are listed below.
Write a python code to solve the above two nonlinear equations, Eqs. (1) and (2), for $E$ with $V_0 = 10$.
(Bonus) Substitute each solution of Eq. (1) $E$ into the symmetric wave functions:
$$
\psi(x)
= \left\{\begin{array}{ll}
C \exp\left(\pi x\sqrt{V_0 - E}\right) + E, & x < -1/2 \\
\cos\left(\pi x\sqrt{E}\right) + E, & -1/2 \le x \le 1/2 \\
C \exp\left(-\pi x\sqrt{V_0 - E}\right) + E, & x > 1/2
\end{array}\right.,
$$
where
$$
C \equiv \frac{\cos\left(\frac{\pi}{2}\sqrt{E}\right)}{\exp\left(-\frac{\pi}{2}\sqrt{V_0-E}\right)}.
$$
Substitute each solution of Eq. (2) $E$ into the asymmetric wave functions:
$$
\psi(x)
= \left\{\begin{array}{ll}
-S \exp\left(\pi x\sqrt{V_0 – E}\right) + E, & x < -1/2 \\
\sin\left(\pi x\sqrt{E}\right) + E, & -1/2 \le x \le 1/2 \\
S \exp\left(-\pi x\sqrt{V_0 - E}\right) + E, & x > 1/2
\end{array}\right.,
$$
where
$$
S \equiv \frac{\sin\left(\frac{\pi}{2}\sqrt{E}\right)}{\exp\left(-\frac{\pi}{2}\sqrt{V_0-E}\right)}.
$$
Plot all of these $\psi(x)$. Note that each solution $E$ to Eq. (1) or (2) leads to a $\psi(x)$ curve.
(15 + (bonus 5) marks)
In [ ]: