lab10¶

`f(x)`¶

Implement a function \(f(x)\) which computes:

\[f(x) = \frac{\exp(-x^2)}{1+x^2} + \frac{2 \cos(x)^2}{1 + (x - 4)^2}\]

Examples

In [ ]: f(0)
Out[ ]: 1.1176470588235294

In [ ]: f(2)
Out[ ]: 0.07293440360502447

While you do not need to do the following, it is possible to implement the function f so that it can compute the function for array input arguments. If you do this, the examples above should work as shown, and in addition you should be able to do:

In [ ]: import numpy as np

In [ ]: x = np.array([0, 2])

In [ ]: f(x)
Out[ ]: array([ 1.11764706,  0.0729344 ])

For the following tasks, it is not necessary to support numpy array arguments x as shown here but you may find this useful to plot \(f(x)\) more conveniently.

Exercise `integrate_f_from0(b)`¶

Implement a function integrate_f_from0(b) which computes (an approximation of)

\[\int\limits_0^bf(x)\mathrm{d} x\]

where the function \(f(x)\) is as defined in the exercise above. The function integrate_f_from0(b) should return a floating point number.

`find_max_f()`¶

Implement a function find_max_f() which returns (an approximation of) the \(x\) for which \(f(x)\) takes the maximum value. Your function find_max_f should return a floating point number.

Exercise `find_f_equals_1()`¶

Implement a function find_f_equals_1() which returns a float which is (an approximation of) the \(x\) for which \(f(x)=1\). We are only interested in the solution where \(x\) is negative.

Exercise `lin_int(xs, ys)`¶

Implement a function lin_int(xs, ys) which accepts sequences xs and ys of data as input, and returns a callable object f(x) which returns y(x) using linear interpolation between the points provided by the data xs and ys.

You can write your own interpolation function or use a library function for this.

Example

In [ ]: xs = [1, 2, 3]

In [ ]: ys = [10, 20, 15]

In [ ]: f = lin_int(xs, ys)

In [ ]: f(1)
Out[ ]: array(10.0)

In [ ]: f(2)
Out[ ]: array(20.0)

In [ ]: f(1.5)
Out[ ]: array(15.0)

In [ ]: f(2.5)
Out[ ]: array(17.5)

`solve_freefall(ts, v0)`¶

`solve_freefall(ts, v0)` (part 1)¶

Implement the function solve_freefall(ts, v0) which solves the following ordinary differential equation for the velocity v as a function of time t:

\[\frac{\mathrm{d}v}{\mathrm{d}t} = g - \frac{Ac}{m}v|v| \qquad \qquad (1)\]

The ordinary differential equation describes the time derivative of the velocity \(v\) (i.e. the acceleration) which is composed of the gravitational acceleration \(g\) acting downwards and a term acting in the opposite direction of the velocity (\(-v|v|Ac/m\)) which originates from (turbulent) friction with the air. (\(|v|\) represents the absolute value of the velocity \(v\). The product \(v|v|\) is proportional to \(v^2\).)

The constants are

the mass \(m=80\,\mathrm{kg}\) of the object,
a surface area \(A=0.5\,\mathrm{m^2}\),
\(c=0.25\,\mathrm{kg/m^3}\) is a friction coefficient, and
\(g=9.81\,\mathrm{m/s^2}\) the Earth’s gravitational acceleration.
This approximates a human-sized object in free fall in the Earth’s atmosphere.

The argument ts is a sequence containing time steps for which we desire the solution \(v(t)\). The first value in ts must be the initial time \(t_0\). The second argument v0=\(v_0\) is a number that represents the initial value for the solution at time \(t_0\), i.e. \(v_0=v(t_0)\).

The function solve_freefall should return an numpy array vs containing (an approximation of) the solution for \(v(t)\) for all elements in ts, i.e. vs=v(ts). For the actual numerical integration use the solve_ivp function from scipy.integrate.

(This exercise is analogous to the way that the function solve_ode0(ts, y0) in the example below computes and returns the numerical solution of the differential equation \(\frac{\mathrm{d}y}{\mathrm{d}t}=-2y\) in the example below.)

Use settings atol=1e-8 and rtol=1e-8 for the numerical integration in the solve_ivp command.

Example:

In [ ]: solve_freefall([0, 1], 0)
Out[ ]: array([0.        , 9.76018245])

In [ ]: solve_freefall([0, 0.5, 1], 0)
Out[ ]: array([0.        , 4.89874415, 9.76018245])

In [ ]: solve_freefall([0, 1], v0=10)
Out[ ]: array([10.        , 19.45769981])

You should plot the solution v(t) to check your implementation — do you understand and can you explain how the velocity changes?

You can also compute the final (asymptotic) velocity analytically when turbulent friction and gravitational pull are equal (in this case dv/dt=0). Does this match the velocity that solve_freefall returns for large times \(t\)?

`solve_freefall(ts, v0, m)` (part 2)¶

Extend or modify your function solve_freefall so that it takes a third argument for the mass \(m\) with default value 80kg.

Example:

In [ ]: solve_freefall([0, 1], v0=0)
Out[ ]: array([0.        , 9.76018245])

In [ ]: solve_freefall([0, 1], v0=0, m=80)
Out[ ]: array([0.        , 9.76018245])

In [ ]: solve_freefall([0, 1], v0=0, m=200)
Out[ ]: array([0.        , 9.78999986])

How does the final velocity change if you change the mass \(m\) of the object? Is this in agreement with the analytical expression for the final velocity that you have used in the previous question?

Appendix: Example solving ODE with `scipy.integrate.solve_ivp`¶

import matplotlib.pyplot as plt
import scipy.integrate
import numpy


def solve_ode0(ts, y0):
    """Solve dy/dt = -2*y with y0=y(ts[0]).

    Returns array ys with (approximation of) solution for all elements in ts, i.e.
    ys=y(ts)
    """
    c = 2  # in dy/dt = -c*x, we have c=2 here

    def rhs(t, y):
        """Compute the right hand side of dy/dt=-c*y."""
        return -c * y

    sol = scipy.integrate.solve_ivp(
        fun=rhs,
        t_span=[min(ts), max(ts)],  # create t_span from ts
        y0=[y0],  # y0 on the right is number, and we need a sequence
        t_eval=ts,  # return t and y for t in ts
        atol=1e-6,  # optional: absolute tolerance
        rtol=1e-6,  # optional: relative tolerance
    )

    return sol.y[0]


def plot_example():
    """Create plot with solution."""
    ts = numpy.linspace(0, 10, 100)
    ys1 = solve_ode0(ts, y0=1)
    plt.plot(ts, ys1, "g-")
    plt.xlabel("t")
    plt.ylabel("y(t)")
    plt.show()


if __name__ == "__main__":
    ts = [0, 0.5, 1]
    ys = solve_ode0(ts, y0=1)
    print("With y0=1, the solution is")
    for t, y in zip(ts, ys):
        print(f"y({t:f})={y:f}")

    plot_example()

Please submit your file lab10.py for this assignment.

Additional (voluntary) tasks are available in lab10-extra.

End of lab10.

lab10¶

f(x)¶

Exercise integrate_f_from0(b)¶

find_max_f()¶

Exercise find_f_equals_1()¶

Exercise lin_int(xs, ys)¶

solve_freefall(ts, v0)¶

solve_freefall(ts, v0) (part 1)¶

solve_freefall(ts, v0, m) (part 2)¶

Appendix: Example solving ODE with scipy.integrate.solve_ivp¶