<jupyter_start><jupyter_text>![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.Held at Technische Universität Bergakademie Freiberg.Comments to:*Prof. Dr. <NAME>  Chair of Soil Mechanics and Foundation Engineering  Geotechnical Institute  Technische Universität Bergakademie Freiberg.*https://tu-freiberg.de/en/soilmechanics  Exercise 3 - Generalized midpoint rule, Picard iterationsThe ODE $\dot{y} = k (y_\text{eq} - y)^2$ is to be integrated numerically.First, we import some libraries for numerical basics and for plotting.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

#Some plot settings
plt.style.use('seaborn-deep')
plt.rcParams['lines.linewidth']= 2.0
plt.rcParams['lines.color']= 'black'
plt.rcParams['legend.frameon']=True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize']=14
plt.rcParams['font.size'] = 14
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
plt.rcParams['figure.figsize'] = (8, 6)
y_eq = 1.
k = 2.
def y_analytical(t):#assuming y0=0
    return k*t*y_eq**2./(1.+k*t*y_eq)<jupyter_output><empty_output><jupyter_text>According to the ODE definition $\dot{y}=f[y(t),t]$ we define the rate function itself:<jupyter_code>def dy_dt(y):
    return k*(y_eq - y)**2<jupyter_output><empty_output><jupyter_text>Let's plot the function and its rate for some exemplary values ($y_0 = 0$) in a time interval $t \in [0,5]$:<jupyter_code>time = np.linspace(0,5,100)
plt.plot(time,y_analytical(time),label='$y$')
plt.plot(time,dy_dt(y_analytical(time)),label='$\\dot{y}$')
plt.legend()
plt.xlabel('$t$ / s')
plt.ylabel('$y$, $\\dot{y}$ / s$^{-1}$')<jupyter_output><empty_output><jupyter_text>Now we want to get the same result with a generalized midpoint rule. It can be considered as a weighted average between an explicit forward Euler and an implicit backward Euler scheme. For $\alpha > 0$ the update proceeds implicitly with the *recursive* update formula$y_{n+1}^{i+1} = y_n + [(1-\alpha) f(y_n) + \alpha f(y_{n+1}^i)] (t_{n+1} - t_n)$We define the corresponding Euler update as$\Delta y_\text{MP}^{i} = y_{n+1}^{i+1} - y_n = [(1-\alpha) f(y_n) + \alpha f(y_{n+1}^i)] \Delta t$<jupyter_code>def Mid_Point_Update(y_i,y_n,al,dt):
    return ((1.-al)*dy_dt(y_n) + al*dy_dt(y_i))*dt<jupyter_output><empty_output><jupyter_text>Again, the first step is a forward Euler step. After that, a mixture of BE and FE schemes occurs.We now need to define the time stepping itself. So we choose a time-step size and the same starting value (initial condition) of $y_0 = y(t=0) = 0$ to initialize our solution field. Now, for every time step we append to this array the value of the next time step by applying the generalized Midpoint update rule. Each new time step value then serves as the initial condition for the following time step.We will *recursively* re-fresh the update with a new estimate for the rate until the estimated solution does not change anymore, i.e.$\left| y_{n+1}^{i+1} - y_{n+1}^{i} \right| < \epsilon_\text{abs} \cup \left| \frac{y_{n+1}^{i+1} - y_{n+1}^{i}}{y_{n+1}^{i}} \right| < \epsilon_\text{rel}$with suitably chosen absolute and/or relative tolerances.For simplicity, let's set only an absolute tolerance of $10^{-6}$ here.<jupyter_code>def integrate(dt,alpha):
    #Startwerte
    t_end = 5.
    absolute_tolerance = 1.e-6
    max_iter = 1000
    iteration_counter = np.array([0])
    y = np.array([0.])
    times = np.array([0.])
    #
    while times[-1]+dt < t_end: #repeat the loop as long as the final time step is below the end point
        times = np.append(times,times[-1]+dt) #here define the next time point as the previous time point plus the time increment dt
        y_old = y[-1] #Starting value for recursive update
        i = 0
        #
        while True:
            y_new = y[-1]+Mid_Point_Update(y_old,y[-1],alpha,dt) #one backward Euler step
            i += 1
            #if change is below tolerance, above maximum iterations, or scheme is fully explicit, stop iterations
            if (np.abs(y_new - y_old) < absolute_tolerance or i > max_iter or alpha == 0.0): 
                break
            y_old = y_new #preparation of next recursion

        y = np.append(y,y_new) #append the new found solution to the solution vector
        iteration_counter = np.append(iteration_counter,i) #store how much iterations this time step took to converge
    return times, y, iteration_counter
from ipywidgets import widgets
from ipywidgets import interact

@interact(dt=widgets.BoundedFloatText(value=0.2,min=1e-3,max=1,step=0.1,description='$\Delta t$ / s'),
          alpha=widgets.RadioButtons(options=[0,0.5,1], value=1, description='$\\alpha$'))

def plot_CN(dt=0.1,alpha=1):
    time_n, y_n, iter = integrate(dt,alpha)
    fig, ax = plt.subplots(ncols=2,nrows=2,figsize=(18,12))
    #solution
    ax[0][0].plot(time,y_analytical(time),label='analytical')
    ax[0][0].plot(time_n,y_n,label='numerical',marker='d',ls=':')
    ax[0][0].legend()
    ax[0][0].set_xlabel('$t$ / s')
    ax[0][0].set_ylabel('$y$')
    #error
    ax[0][1].plot(time_n,y_n-y_analytical(time_n),ls=':',marker='d',label='absolute')
    ax[0][1].plot(time_n[1:],(y_n[1:]-y_analytical(time_n[1:]))/y_analytical(time_n[1:]),
               ls=':',marker='d',label='relative')
    ax[0][1].legend()
    ax[0][1].set_xlabel('$t$ / s')
    ax[0][1].set_ylabel('error')
    #iterations
    ax[1][0].plot(iter)
    ax[1][0].set_xlabel('time step')
    ax[1][0].set_ylabel('number of iterations')
    
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15280165005156413, 'fraction_numerical': 0.021656926778961843, 'mean_word_length': 4.48444863336475, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 24, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '324910', 'license': 'MIT', 'n_tokens_mistral': 2167, 'n_tokens_neox': 2052, 'n_words': 524, 'path': '03b_Generalized_Midpoint_Picard_Nonlin.ipynb', 'repo_name': 'dominik-kern/Numerical_Methods_Introduction', 'size': '5824'}