Chapter 5

Plot Taylor polynomial approximations to \(\sin x\). The sine function can be approximated by a polynomial according to the following formula: $$ \sin x \approx S(x ; n)=\sum_{j=0}^{n}(-1)^{j} \frac{x^{2 j+1}}{(2 j+1) !} $$ The expression \((2 j+1) !\) is the factorial (see Exercise 3.14). The error in the approximation \(S(x ; n)\) decreases as \(n\) increases and in the limit we have that \(\lim _{n \rightarrow \infty} S(x ; n)=\sin x\). The purpose of this exercise is to visualize the quality of various approximations \(S(x ; n)\) as \(n\) increases. The first part of the exercise is to write a Python function \(\mathrm{S}(\mathrm{x}\), (n) that computes \(S(x ; n)\). Use a straightforward approach where you compute each term as it stands in the formula, i.e., \((-1)^{j} x^{2 j+1}\) divided by the factorial \((2 j+1) !\). (We remark that Exercise A.16 outlines a much more efficient computation of the terms in the series.) The next part of the exercise is to plot \(\sin x\) on \([0,4 \pi]\) together with the approximations \(S(x ; 1), S(x ; 2), S(x ; 3), S(x ; 6)\), and \(S(x ; 12)\) Name of program file: plot_Taylor_sin.py.

Animate a planet's orbit. A planet's orbit around a star has the shape of an ellipse. The purpose of this exercise is to make an animation of the movement along the orbit. One should see a small disk, representing the planet, moving along an elliptic curve. An evolving solid line shows the development of the planet's orbit as the planet moves. The points \((x, y)\) along the ellipse are given by the expressions $$ x=a \cos (\omega t), \quad y=b \sin (\omega t) $$ where \(a\) is the semimajor axis of the ellipse, \(b\) is the semiminor axis, \(\omega\) is an angular velocity of the planet around the star, and \(t\) denotes time. One complete orbit corresponds to \(t \in[0,2 \pi / \omega] .\) Let us discretize time into time points \(t_{k}=k \Delta t\), where \(\Delta t=2 \pi /(\omega n) .\) Each frame in the movie corresponds to \((x, y)\) points along the curve with \(t\) values \(t_{0}, t_{1}, \ldots, t_{i}, i\) representing the frame number \((i=1, \ldots, n)\). Let the plot title of each frame display the planet's instantaneous velocity magnitude. This magnitude is the length of the velocity vector $$ \left(\frac{d x}{d t}, \frac{d y}{d t}\right)=(-\omega a \sin (\omega t), \omega b \cos (\omega t)) $$ which becomes \(\omega \sqrt{a^{2} \sin ^{2}(\omega t)+b^{2} \cos ^{2}(\omega t)}\) Implement the visualization of the planet's orbit using the method above. Run the special case of a circle and verify that the magnitude of the velocity remains constant as the planet moves. Name of program file: planet_orbit.py.

Plot the velocity profile for pipeflow. A fluid that flows through a (very long) pipe has zero velocity on the pipe wall and a maximum velocity along the centerline of the pipe. The velocity \(v\) varies through the pipe cross section according to the following formula: $$ v(r)=\left(\frac{\beta}{2 \mu_{0}}\right)^{1 / n} \frac{n}{n+1}\left(R^{1+1 / n}-r^{1+1 / n}\right) $$ where \(R\) is the radius of the pipe, \(\beta\) is the pressure gradient (the force that drives the flow through the pipe), \(\mu_{0}\) is a viscosity coefficient (small for air, larger for water and even larger for toothpaste), \(n\) is a real number reflecting the viscous properties of the fluid ( \(n=1\) for water and air, \(n<1\) for many modern plastic materials), and \(r\) is a radial coordinate that measures the distance from the centerline \((r=0\) is the centerline, \(r=R\) is the pipe wall). Make a function that evaluates \(v(r) .\) Plot \(v(r)\) as a function of \(r \in[0, R]\), with \(R=1, \beta=0.02, \mu=0.02\), and \(n=0.1 .\) Thereafter, make an animation of how the \(v(r)\) curves varies as \(n\) goes from 1 and down to \(0.01\). Because the maximum value of \(v(r)\) decreases rapidly as \(n\) decreases, each curve can be normalized by its \(v(0)\) value such that the maximum value is always unity. Name of program file: plot_velocity_pipeflow.py.

Plot a w-like function. Define mathematically a function that looks like the 'w' character. Plot this function. Name of program file: plot_w.py.

Plot a smoothed "hat" function. The "hat" function \(N(x)\) defined by (3.5) on page 109 has a discontinuity in the derivative at \(x=1\). Suppose we want to "round" this function such that it looks smooth around \(x=1\). To this end, replace the straight lines in the vicinity of \(x=1\) by a (small) cubic curve $$ y=a(x-1)^{3}+b(x-1)^{2}+c(x-1)+d $$ for \(x \in[1-\epsilon, 1+\epsilon]\), where \(a, b, c\), and \(d\) are parameters that must be adjusted in order for the cubic curve to match the value and the derivative of the function \(N(x)\). The new rounded functions has the specification $$ \tilde{N}(x)= \begin{cases}0, & x<0 \\ x, & 0 \leq x<1-\epsilon \\\ a_{1}(x-1)^{3}+b(x-1)+c(x-1)+d_{1}, & 1-\epsilon \leq x<1 \\\ a_{2}(x-1)^{3}+b(x-1)+c(x-1)+d_{2}, & 1 \leq x<1+\epsilon \\ 2-x, & 1+\epsilon \leq x<2 \\ 0, & x \geq 2\end{cases} $$ with \(a_{1}=\frac{1}{3} \epsilon^{-2}, a_{2}=-a_{1}, d_{1}=1-\epsilon+a_{1} \epsilon^{3}, d_{2}=1-\epsilon-a_{2} \epsilon^{3}\), and \(b=c=0 .\) Plot this function. (Hint: Be careful with the choice of \(x\) coordinates!) Name of program file: plot_hat.py.

Experience overflow in a function. When an object (ball, car, airplane) moves through the air, there is a very, very thin layer of air close to the object's surface where the air velocity varies dramatically \(^{18}\), from the same value as the velocity of the object at the object's surface to zero a few centimeters away. The change in velocity is quite abrupt and can be modeled by the functiion $$ v(x)=\frac{1-e^{x / \mu}}{1-e^{1 / \mu}} $$ where \(x=1\) is the object's surface, and \(x=0\) is some distance away where one cannot notice any wind velocity \(v\) because of the passing object \((v=0)\). The vind velocity coincides with the velocity of the object at \(x=1\), here set to \(v=1\). The parameter \(\mu\) is very small and related to the viscosity of air. With a small value of \(\mu\), it becomes difficult to calculate \(v(x)\) on a computer. Make a function \(v(x, m u=1 E-6\), exp=math.exp) for calculating the formula for \(v(x)\) using exp as a possibly user-given exponentional function. Let the v function return the nominator and denominator in the formula as well as the fraction (result). Call the v function for various \(\mathrm{x}\) values between 0 and 1 in a for loop, let mu be \(1 \mathrm{E}-3\), and have an inner for loop over two different exp functions: math. exp and numpy. exp. The output will demonstrate how the denominator is subject to overflow and how difficult it is to calculate this function on a computer. Also plot \(v(x)\) for \(\mu=1,0.01,0.001\) on \([0,1]\) using 10,000 points to see what the function looks like. Name of program file: boundary_layer_func1.py.