Innlevering 3

Du leverer besvarelsen din som to filer:

Oppgavene 1-3: Finnes i Oppgave filen på wikisiden og besvarelsen leveres som en pdf fil gjennom Ovsys2.

Oppgave 4: Finnes her. Legg til svarene dine og lever Jupyter notebook gjennom Ovsys2.

4) The goal of this exercise is to approximate solutions of two nonlinear systems of ODEs from the book:

-) The population dynamics: \begin{align*} \frac{d y_1}{dt}& = C_1y_1 (1-b_1y_1-d_2y_2),\\ \frac{dy_2}{dt}&= -C_2y_2(1-b_2y_2-d_1y_1), \end{align*} where $y_1$ and $y_2$ denote the two competing populations, $C_1$ and $C_2$ represent the growth rates of the two populations, $d_1$ and $d_2$ govern the interaction between the two populations and $b_1$ and $b_2$ are related to the available quantity of nutrients.

-) The Baseball trajectory: \begin{align*} \frac{d{\bf x}}{dt}& = {\bf v},\\ \frac{d{\bf v}}{dt}& = {\bf F}, \end{align*} where ${\bf x}(t)=(x(t), y(t), z(t))^{T}$ denotes the position of the ball at time $t$ and ${\bf v}(t)= (v_x(t), v_y(t), v_z(t))^{T}$ its velocity. Moreover, the vector ${\bf F}$ is given by \begin{align*} F_x& =-F(v)vv_x+B\omega(v_z\sin(\phi)-v_y\cos(\phi)),\\ F_y& =-F(v)vv_y+B\omega v_x\cos(\phi),\\ F_z&= -g-F(v)vv_z-B\omega v_x\sin(\phi), \end{align*} where \begin{equation*} F(v)= 0.0039+\frac{0.0058}{1+e^{\frac{v-35}{5}}} \end{equation*} is a friction coefficient, $g=9.8$, $v$ denotes the modulus of ${\bf v}$, $B= 4.1 \cdot10^{-4}$, $\phi=1.4$ degree is the pitching angle, and $\omega=180\cdot 1,047198$ degree per second, the modulus of the angular velocity impressed to the ball from the pitcher.

a) Read the documentation of solve_ivp to understand exactly how to implement the above ODEs in the format required for the use of the function "solve_ivp" of scipy using the option "method='RK45'.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

Here we implemented the population dynamics problem. Use this as an exampe for implementing the first order ODE system describing the Baseball projectory.

#ODE problems

# Population dynamics
C1=1
C2=1
d1=1
d2=1
b1=0
b2=0
def populationdyn(t, y):
    y0, y1 = y
    f0 = C1 * y0 * (1 - b1 * y0 - d2 * y1)
    f1 = -C2 * y1 * (1 - b2 * y1 - d1 * y0)
    return np.array([f0, f1])

# Baseball trajectory 
def baseball(t, y):
    # Write your code here
    return

b) The function "generateA" provided her can be used to generate a random $n\times n$ matrix $A$, which is guaranteed to be diagonalizable. This allows us to define a linear system of ODEs $y'=Ay$, $y(t_0)=y_0$.

#ODE problems

# Linear system of ODEs
# We generate randomly a linear system of ODEs y'=Ay, where we make sure 
# that A is a diagonalizable matrix
# Generate a random diagonalizable matrix
def generateA(n):
    rank = 0
    while rank < n:
        A = -np.random.rand(n,n)
        # the following  function of numpy.linalg computes eigenvalues and eigenvectors 
        eigenvaluesA, eigenvectorsA = np.linalg.eig(A)
        # The following function of numpy.linalg estimates the rank
        rank = np.linalg.matrix_rank(eigenvectorsA)
        print("Estimated rank of the matrix of egenvectors:", rank)
    return A, eigenvaluesA, eigenvectorsA

# We next use the matrix to define the ODE y'=Ay
def linearODE(t,y):
    #y vector with n components
    f = A@y
    return f

Examine the above code and explain briefly what the code of "generateA" does. In particular, explain why the resulting $n\times n$ matrix $A$ is guaranteed to be diagonalizable.

Write your answer here.

c) Use the code "solve_ivp" of scipy.integrate with relative tolerance rtol=2.3*1e-14, and absolute tolerance atol=1e-15 to solve the following 3 ODE problems:

1. the population dynamics system on the time interval $[0,10]$ with initial value y0= np.array([2.0 , 2.0]). Furthermore, make a phase-space plot, i.e., plot the solution in the plane with y[0] on the x-axis and y[1] on the y-axis. What do you observe?
2. the Baseball trajectory system on the time interval $[0,10]$ with initial velocity ${\bf v}(0)=v_0(\cos(\phi), 0 \sin(\phi))^{T}$, $v_0=38m/s$. If $\textbf{x}(0)=\textbf{0}$, after how many seconds (approximately) will the ball touch the ground, i.e., $z(t)=0$?
3. the linear ODE system with A generated by "generateA" with $n=10$, $20$ and $30$, on the time interval $[0,1]$.

For each system plot the numerical solution obtained by "solve_ivp" on the given time intervals.

# Approximations to the population dynamics equations with solve_ivp
y0= np.array([2.0 , 2.0])

T=10.
sol = solve_ivp(populationdyn, [0, T], y0, method='RK45', rtol=2.3*1e-14, atol=1e-15)
plt.plot(sol.t,sol.y[0],'k',
        sol.t,sol.y[1],'b')

#Phase space plot: we plot the solution in the plane with y[0] on the x-asis and y[1] on the y-asis
# Write your plot command here

Your answer to: What did you observe?

# Approximation to the baseball trajectory with solve_ivp 
# Write your code here

# Write your plot command here

Your answer to: After how many seconds does the ball touch the ground?

# Linear ODE. we use also here solve_ivp and plot the solution. We use a random initial vector with n components
n=100
A, eigenvaluesA, eigenvectorsA = generateA(n)
y0 = np.random.rand(n)
T=1.
# Solve the linearODE IVP with "solve_ivp" here with y0 as initial vector and T=1. as final time
sol = solve_ivp(linearODE, [0, T], y0, method='RK45', rtol=2.3*1e-14, atol=1e-15)
# Then plot the first, second and 10th component of the obtained numerical solution here

How accurate do you expect the solution to be when using tolerances rtol=2.3*1e-14, and atol=1e-15? What does it mean that you use " method='RK45' "? Read the documentation of "solve_ivp" and provide (short) answers.

Write your answer here.

d) The function "generateA" also provides approximations of the eigenvalues and eigenvectors of the matrix $A$ using the routine "numpy.linalg.eig", which computes accurate approximations of eigenvalues and eigenvectors of a given matrix $A$ of moderate size.

Use this information to provide a code, which uses the method described in Section 6.9 in the book A. Hole, Kalkulus and linear algebra, Bind 2, that generates approximations of the solution of the linear ODE problem at times $[t_0,t_1,\dots ,t_N]$ where $t_N=T$ is the final time of the interval of integration $[0,T]$, $t_{k+1}-t_k=h$, and $h=T/N$, where $N$ is the total number of time values. For comparison, plot in the same plot the numerical solution you obtain with this method and the numerical solution you obtained earlier with "solve_ivp".

Furthermore, compare also the two solutions at final time $T$ by computing the norm of the difference of the two numerical solutions. Make sure to use the same initial vector and print the value of the norm of the difference.

# Now we solve the same linear ODE with the method that we learned in class using the factorization of A = U D U^{-1} 
# where D is a diagonal matrix with the eigenvalues of A along the diagonal and U the matrix whose columns are the 
# corresponding eigenvectors
# the egenvectors of A are in eigenvectorsA

# N = total number of time values 
N=100;
E=np.zeros((y0.size,N+1))
# sollinODE is the array where you will store the values of the numerical solution at times t0=0,t_1,t_2, ...t_N=T, 
# an array with the same number of rows as y0 and number of columns = N+1
sollinODE=np.zeros((y0.size,N+1),dtype=complex)
# time values t0=0,t_1,t_2, ...t_N=T with t_n-t_{n-1}=T/N
tt = np.linspace(0, T, N+1)
# Let U=eigenvectorsA
# D=np.diag(eigenvaluesA)
# At time t the solution is now U exp(tD) U^{-1}y0
# exp(tD) is the diagonal matrix with exp(t\lambda_k) k=1:n along the diagonal, and \lambda_k is an eigenvalue of A
# Use this formula to write the solution for t0=0,t_1,t_2, ...t_N=T,
#
sollinODE[:,0] = y0
for n in range(N):
    # put here the lines of code needed for implementing the method



#Compare with the numerical solution obtained with solve_ivp in a plot. Use marker "+" for the sollinODE solution 
# and marker "-" for the soluition given by "solv_ivp". Use different colors for the different components of the solution.
# Plot three components of your choice
#
# Insert the plot command for the solution solllinODE here 

plt.plot(sol.t,sol.y[0],'k-',
        sol.t,sol.y[1],'b-',
        sol.t,sol.y[9],'r-')
#Compare with the numerical solution obtained with solve_ivp: compute the norm of the difference of the two solution at final time T
# Use the following two lines of code:
norm_diff = np.linalg.norm(sollinODE[:,-1]-sol.y[:,-1])
print("Norm of difference between the two numerical approximations at time T:", norm_diff)

e) Finally, you will have to implement your own version of the following method of order $3$ invented by Kutta: \begin{align*} y_{n+1}&=y_n+\frac{h}{6}(k_1+4k_2+k_3)\\ k_1& =f(t_n,y_n)\\ k_2& =f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)\\ k_3& =f(t_n+h,y_n-hk_1+2hk_2). \end{align*}

You find below an implementation of the Forward Euler and of the Improved Euler method. Make your implementation of the Kutta method in the same format as used for these two methods.

# We now implement the forward Euler method (order 1), the Improved Euler method (order 2) and the Kutta method (order 3) for a system of odes as defined in cell 1
# 
def FE(odefile, iv, h, N):
    solFE=np.zeros((iv.size,N+1))
    solFE[:,0]= iv
    t0=0.0
    tt=np.zeros(N+1)
    for n in range(N):
        tn=t0+n*h   
        sol_o=solFE[:,n]
        k1=odefile(tn,sol_o)
        sol_n=sol_o+h*k1
        solFE[:,n+1]=sol_n
        tt[n+1]=tn+h
    return solFE, tt
#    
def IE(odefile, iv, h, N):
    solIE=np.zeros((iv.size,N+1))
    solIE[:,0]= iv
    t0=0.0
    tt=np.zeros(N+1)
    for n in range(N):
        tn=t0+n*h   
        sol_o=solIE[:,n]
        k1=odefile(tn,sol_o)
        k2=odefile(tn+h/2,sol_o+h/2*k1)
        sol_n=sol_o+h*k2
        solIE[:,n+1]=sol_n
        tt[n+1]=tn+h
    return solIE, tt
#
#    
def kutta(odefile, iv, h, N):
    sol_kutta=np.zeros((iv.size,N+1))
    sol_kutta[:,0]= iv
    t0=0.0
    tt=np.zeros(N+1)
    for n in range(N):
        tn=t0+n*h 
        #insert you code here ...
        
    return sol_kutta, tt
#

Next, design a numerical experiment that checks numerically the order of the methods for each of the three ODE problems: the population dynamics, the baseball trajectory and linear ODE systems.

• Run multiple experiments on a fixed time interval $[0,T]$ (take for example $T=1$) with several different step-sizes (e.g. $h=1/2,1/4,1/8,1/16,\dots$).
• Record the errors given by the methods at time $T$. Therefore you will have to use "solve_ivp" to obtain an accurate reference solution. Use low tolerances, e.g. rtol=2.3*1e-14, and atol=1e-15.
• Then plot the errors you obtained for the various methods as a function of $h$ in a loglog-plot using the command plt.loglog, which uses logarithmic scale both on the x-axis and on the y-axis. If the methods are implemented correctly, you will see that in a loglog-plot the errors decrease linearly with slope $1$ for the Forward Euler (order $1$), $2$ for the Improved Euler (which has order $2$) and $3$ for the Kutta method (of order $3$). SInce the order of the methods is $\mathcal{O}(h^p)$, you need to use "small enough" values of $h$ to observe the correct order.

Do you observe order $3$ for the Kutta method?

#Order test experiment you can use oderfile = linearODE or baseball or populationdyn
odefile = populationdyn
# provide initial value of correct dimension (compatible with odefile)
y0= np.array([1.0 , 0.2 ])
#y0 = np.random.rand(100)
# y0 = 
# final time
T=1.
# Compute reference solution with "solve_ivp" with low rel and abs tolerance
sol = solve_ivp(odefile, [0, T], y0, method='RK45', rtol=1e-14, atol=1e-15)
maxiter= 100
K=8
eFE=np.zeros(K)
eIE=np.zeros(K)
ekutta=np.zeros(K)
hh=np.zeros(K)
for k in range(K):
    N=2**(k+1)
    h=T/N
    hh[k]=h
    solFE = np.zeros((y0.size,N+1))
    solIE = np.zeros((y0.size,N+1))
    sol_kutta =  np.zeros((y0.size,N+1))
    solFE, tt = FE(odefile, y0, h, N)
    eFE[k]=np.linalg.norm(sol.y[:,-1]-solFE[:,-1])
    solIE, tt = IE(odefile, y0, h, N)
    eIE[k]=np.linalg.norm(sol.y[:,-1]-solIE[:,-1])
    # Add here the lines of code to check the order of the method of Kutta
    sol_kutta, tt = kutta(odefile, y0, h, N)
    ekutta[k]=np.linalg.norm(sol.y[:,-1]-sol_kutta[:,-1])

plt.loglog(hh,eFE,'b-o',
          hh,eIE,'r-o',
          hh,ekutta,'g-o',           
          hh,hh**(1),'k-o',
          hh,hh**(2),'k-o',
          hh,hh**(3),'k-o' )    

Your answer to: Do you observe order 3 for the Kutta method?