Bone remodeling

Problem setting

Import modules :

import numpy as np import matplotlib.pyplot as plt

Coefficients :

## DO NOT CHANGE THESE : a1= 3 a2= 4 b1= 0.2 b2= 0.02 g11= 1.1 g12= 1 g21=-0.5 g22= 0

Initial condition :

# initial condition : DO NOT MODIFY THIS y0 = [1.060660171777250 +10, 2.121320343560527e+02]

Time stepping options :

## ONLY CHANGE Ntstep Ntsteps = 100000 t = np.linspace(0, 1000, num=Ntsteps+1) # DO NOT CHANGE THIS PART

Part 1 : Solving the problem

Implement and solve the problem using your own functions for the backward Euler's method and Newton's method.

import numpy as np def Newton(y_0,g,J,Nmax=1,eps1=1e-1,eps2=1e-1,disp=True): y = y_0 Delta = [] R = [] gk = g(y) Jk = J(y) r = np.linalg.norm(gk) dy = np.linalg.solve(Jk,gk) delta = np.linalg.norm(dy) R.append(r) Delta.append(delta) k = 0 if disp: print('residual : ',r,' -- Increment : ',delta) while (k<Nmax) and (delta > eps1) and (r > eps2): k+=1 y -= dy gk = g(y) Jk = J(y) r = np.linalg.norm(gk) dy = np.linalg.solve(Jk,gk) delta = np.linalg.norm(dy) R.append(r) Delta.append(delta) if disp: print('After',k,'iteration(s) = residual : ',r,' -- Increment : ',delta) # write code here return y-dy, Delta, R

def f(t,y): return np.array([a1*(y[0]**g11)*(y[1]**g21)-(b1*y[0]), a2*(y[0]**g12)*(y[1]**g22)-(b2*y[1])]) def L(t,y): return np.array([[a1*g11*(y[0]**(g11-1))*(y[1]**g21)-b1, g21*a1*(y[0]**g11)*(y[1]**(g21-1))],[g12*a2*(y[0]**(g12-1)*(y[1]**g22)), g22*a2*(y[0]**g12)*(y[1]**(g22-1))-b2]]) def BackwardEuler_M(tk,yk,f,L,M,dt): def g(y): return M.dot(y - yk) - dt*f(tk, y) def J(y): return M-dt*L(tk,y) yk1 = Newton(yk,g,J,disp=False) return yk1[0] M = np.array([[1,0],[0,1]]) y1 = np.zeros(Ntsteps+1) y2 = np.zeros(Ntsteps+1) # Set the initial condition for T y1[0] = y0[0] y2[0] = y0[1] for k in range(Ntsteps): dt = t[k+1]-t[k] yk = np.array([y1[k],y2[k]]) yk1 = BackwardEuler_M(t[k],yk,f,L,M,dt) y1[k+1]=yk1[0] y2[k+1]=yk1[1]

plt.style.use('seaborn-poster') #%matplotlib inline fig, ax1 = plt.subplots() color = 'tab:red' ax1.set_xlabel('Time') ax1.set_ylabel('y1', color=color) plt.ylim(-1,13.2) ax1.plot(t,y1, color) ax1.tick_params(axis='y', labelcolor=color) ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis color = 'tab:blue' ax2.set_ylabel('y2', color=color) # we already handled the x-label with ax1 plt.ylim(-55,850) ax2.plot(t,y2, color) ax2.tick_params(axis='y', labelcolor=color) fig.tight_layout() # otherwise the right y-label is slightly clipped plt.show()

/tmp/ipykernel_160/2476182983.py:1: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-poster')

Part 2 : Convergence Assessment

Compute successive numerical solutions as you increase the number of time steps Ntsteps

You can use 1D piecewise linear interpolation directly from SciPy:

from scipy.interpolate import interp1d

Compute and plot (in double logarithmic scale) the rmse(m) by increasing m until you consider that your solution is accurate enough.

Part 3

Solve the same problem using SciPy's function solve_ivp from the integrate submodule

from scipy.integrate import solve_ivp

Use the integration method method = 'LSODA' . Check what the method does and how to use it.

The methods does not require a time stepping vector. Instead the user sets the absolute tolerance through the keyword argument atol, and the time step size is adapted automatically to guarantee the user-specified accuracy requirements. Compare the number of time steps that are needed by LSODA and your Backward Euler Method to achieve comparable accuracy.

ALEMO_HW2 - Duplicate

Bone remodeling

Problem setting

Part 1 : Solving the problem

Part 2 : Convergence Assessment

Part 3

Appendix

.css-ytumd6{-webkit-text-decoration:none;text-decoration:none;}Bone remodeling

Problem setting

Part 1 : Solving the problem

Part 2 : Convergence Assessment

Part 3

Appendix

Bone remodeling