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()
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.