Report -Bone remodeling
ALEMO HOMEWORK 2
Author: Joël Danton
24.01.2023
Background/Motivations
Problem Statement
Methods
Part 1: Solving the problem with Backwar-Euler's and Newton's method
#################################
# INITIALISATION #
#################################
import numpy as np
import matplotlib.pyplot as plt
#COEFFICIENT ###
## 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
y0 = np.array([1.060660171777250 +10, 2.121320343560527e+02])
### Time stepping options ###
## ONLY CHANGE Ntstep
Ntsteps = 10000
t = np.linspace(0, 1000, num=Ntsteps+1) # DO NOT CHANGE THIS PART
dt=(t[-1]-t[0])/(Ntsteps+1)
y1 = np.zeros(Ntsteps+1)
y2 = np.zeros(Ntsteps+1)
# Set the initial condition
y1[0] = y0[0]
y2[0]= y0[1]
#################################
# NEWTON METHOD #
#################################
def NewtonM(g,J,z0,Nmax=100,eps1=10e-6,eps2=10e-6):
delta_k = np.linalg.norm(np.linalg.solve(J(z0,M,L,dt),g(z0,M,f,dt)))
Res_k= np.linalg.norm(g(z0,M,f,dt))
k=0
zk=z0
zk1=zk-np.linalg.solve(J(z0,M,L,dt),g(z0,M,f,dt))
zk=zk1
while k<=Nmax and delta_k > eps1 and Res_k > eps2:
zk1 -= np.linalg.solve(J(zk,M,L,dt),g(zk,M,f,dt))
delta_k = np.linalg.norm(np.linalg.solve(J(zk,M,L,dt),g(zk,M,f,dt)))
Res_k= np.linalg.norm(g(zk,M,f,dt))
k+=1
zk=zk1
z=zk
return z
#################################
# BACKWARD EULER METHOD #
#################################
M = np.array([[1.,0],[0., 1.]])
def f(y):
Result=[]
Result.append(a1*y[0]**g11*y[1]**g21-b1*y[0])
Result.append(a2*y[0]**g12*y[1]**g22-b2*y[1])
return np.array(Result)
def L(y):
df1_dy1=g11*a1*y[0]**(g11-1)*y[1]**g21-b1
df1_dy2=g21*a1*y[0]**g11*y[1]**(g21-1)
df2_dy1=g12*a2*y[1]**g22
df2_dy2=g22*a2*y[0]**g12*y[1]**(g22-1)-b2
return np.array([[df1_dy1,df1_dy2],[df2_dy1,df2_dy2]])
def BEM(yk,M,f,L,dt):
def g(z,M,f,dt):
return M.dot(z-yk)-dt*f(z)
def J(z,M,L,dt):
return M-dt*L(z)
yk1=NewtonM(g,J,yk)
return yk1
for k in range(0,Ntsteps):
yk = np.array([y1[k],y2[k]])
yk1 = BEM(yk,f,L,M,dt)
y1[k+1]=yk1[0]
y2[k+1]=yk1[1]
Part 2: Convergence assessment
#################################
# LOOP PART 2 #
#################################
Y1=[]
Y2=[]
T=[]
RMSE_Y_m=[]
m_vect=[]
index=0
for j in range (-5000,6000,2000):
m=[Ntsteps+j,2*(j+Ntsteps)]
m_vect.append(int(m[0]))
for i in range (0,2):
t = np.linspace(0, 1000, num=m[i]+1) # DO NOT CHANGE THIS PART
dt=(t[-1]-t[0])/(m[i]+1)
y1 = np.zeros(m[i]+1)
y2 = np.zeros(m[i]+1)
# Set the initial condition
y1[0] = y0[0]
y2[0]= y0[1]
for k in range(0,m[i]):
yk = np.array([y1[k],y2[k]])
yk1 = BEM(yk,f,L,M,dt)
y1[k+1]=yk1[0]
y2[k+1]=yk1[1]
T.append(t) #Save the time stepping
Y1.append(y1) #Save the Y1 values for m and 2m
Y2.append(y2) #Save the Y2 values for m and 2m
#################################
# RMSE PART 2 #
#################################
Y1_m=interpolate.interp1d(T[index],Y1[index], kind='linear')
Y1_m_Size2m=Y1_m(T[index+1])
Y2_m=interpolate.interp1d(T[index],Y2[index], kind='linear')
Y2_m_Size2m=Y2_m(T[index+1])
Y_m=np.array([[Y1_m_Size2m],[Y2_m_Size2m]])
Y_2m=np.array([[Y1[index+1]],[Y2[index+1]]])
rmse_Y_m=np.linalg.norm(Y_m-Y_2m)/(np.sqrt(2*m[0]))
RMSE_Y_m.append(float(rmse_Y_m))
index+=2
Part 3: Solving the problem with SciPy's function
#################################
# Solve with SciPy's function #
#################################
def system(t,state,a,g,b):
y1,y2=state
dy1_dt=a1*y1**g11*y2**g21-b1*y1
dy2_dt=a2*y1**g12*y2**g22-b2*y2
return [dy1_dt,dy2_dt]
p=([a1,a2],[[g11,g12],[g21,g22]],[b1,b2])
time=[0,1000]
res_ex=solve_ivp(system,time,y0,method="LSODA",args=p,atol=0.000001)
Results
Plot part 2: Root mean squared error
Conclusions
Bibliography
Appendix
Code part 1
import numpy as np
import matplotlib.pyplot as plt
#################################
# INITIALISATION #
#################################
#COEFFICIENT ###
## 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
y0 = np.array([1.060660171777250 +10, 2.121320343560527e+02])
### Time stepping options ###
## ONLY CHANGE Ntstep
Ntsteps = 10000
t = np.linspace(0, 1000, num=Ntsteps+1) # DO NOT CHANGE THIS PART
dt=(t[-1]-t[0])/(Ntsteps+1)
y1 = np.zeros(Ntsteps+1)
y2 = np.zeros(Ntsteps+1)
# Set the initial condition
y1[0] = y0[0]
y2[0]= y0[1]
#################################
# NEWTON METHOD #
#################################
def NewtonM(g,J,z0,Nmax=100,eps1=10e-6,eps2=10e-6):
delta_k = np.linalg.norm(np.linalg.solve(J(z0,M,L,dt),g(z0,M,f,dt)))
Res_k= np.linalg.norm(g(z0,M,f,dt))
k=0
zk=z0
zk1=zk-np.linalg.solve(J(z0,M,L,dt),g(z0,M,f,dt))
zk=zk1
while k<=Nmax and delta_k > eps1 and Res_k > eps2:
zk1 -= np.linalg.solve(J(zk,M,L,dt),g(zk,M,f,dt))
delta_k = np.linalg.norm(np.linalg.solve(J(zk,M,L,dt),g(zk,M,f,dt)))
Res_k= np.linalg.norm(g(zk,M,f,dt))
k+=1
zk=zk1
z=zk
return z
#################################
# BACKWARD EULER METHOD #
#################################
M = np.array([[1.,0],[0., 1.]])
def f(y):
Result=[]
Result.append(a1*y[0]**g11*y[1]**g21-b1*y[0])
Result.append(a2*y[0]**g12*y[1]**g22-b2*y[1])
return np.array(Result)
def L(y):
df1_dy1=g11*a1*y[0]**(g11-1)*y[1]**g21-b1
df1_dy2=g21*a1*y[0]**g11*y[1]**(g21-1)
df2_dy1=g12*a2*y[1]**g22
df2_dy2=g22*a2*y[0]**g12*y[1]**(g22-1)-b2
return np.array([[df1_dy1,df1_dy2],[df2_dy1,df2_dy2]])
def BEM(yk,M,f,L,dt):
def g(z,M,f,dt):
return M.dot(z-yk)-dt*f(z)
def J(z,M,L,dt):
return M-dt*L(z)
yk1=NewtonM(g,J,yk)
return yk1
for k in range(0,Ntsteps):
yk = np.array([y1[k],y2[k]])
yk1 = BEM(yk,f,L,M,dt)
y1[k+1]=yk1[0]
y2[k+1]=yk1[1]
#################################
# PLOT PART 1 #
#################################
fig, ax1 = plt.subplots()
plt.title("Evolution in time of $y_1$ and $y_2$\nPART 1")
color = 'tab:red'
ax1.set_xlabel('Time $t$ [s]')
ax1.set_ylabel('$y_1$', color=color)
ax1.plot(t,y1, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()
color = 'tab:blue'
ax2.set_ylabel('$y_2$', color=color)
ax2.plot(t,y2,color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()
plt.show()
Code part 2
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
#################################
# INITIALISATION #
#################################
#COEFFICIENT ###
## 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
y0 = np.array([1.060660171777250 +10, 2.121320343560527e+02])
### Time stepping options ###
## ONLY CHANGE Ntstep
Ntsteps = 10000
#################################
# NEWTON METHOD #
#################################
def NewtonM(g,J,z0,Nmax=100,eps1=10e-6,eps2=10e-6):
delta_k = np.linalg.norm(np.linalg.solve(J(z0,M,L,dt),g(z0,M,f,dt)))
Res_k= np.linalg.norm(g(z0,M,f,dt))
k=0
zk=z0
zk1=zk-np.linalg.solve(J(z0,M,L,dt),g(z0,M,f,dt))
zk=zk1
while k<=Nmax and delta_k > eps1 and Res_k > eps2:
zk1 -= np.linalg.solve(J(zk,M,L,dt),g(zk,M,f,dt))
delta_k = np.linalg.norm(np.linalg.solve(J(zk,M,L,dt),g(zk,M,f,dt)))
Res_k= np.linalg.norm(g(zk,M,f,dt))
k+=1
zk=zk1
z=zk
return z
#################################
# BACKWARD EULER METHOD #
#################################
M = np.array([[1.,0],[0., 1.]])
def f(y):
Result=[]
Result.append(a1*y[0]**g11*y[1]**g21-b1*y[0])
Result.append(a2*y[0]**g12*y[1]**g22-b2*y[1])
return np.array(Result)
def L(y):
df1_dy1=g11*a1*y[0]**(g11-1)*y[1]**g21-b1
df1_dy2=g21*a1*y[0]**g11*y[1]**(g21-1)
df2_dy1=g12*a2*y[1]**g22
df2_dy2=g22*a2*y[0]**g12*y[1]**(g22-1)-b2
return np.array([[df1_dy1,df1_dy2],[df2_dy1,df2_dy2]])
def BEM(yk,M,f,L,dt):
def g(z,M,f,dt):
return M.dot(z-yk)-dt*f(z)
def J(z,M,L,dt):
return M-dt*L(z)
yk1=NewtonM(g,J,yk)
return yk1
#################################
# LOOP PART 2 #
#################################
Y1=[]
Y2=[]
T=[]
RMSE_Y_m=[]
m_vect=[]
index=0
for j in range (-5000,6000,2000):
m=[Ntsteps+j,2*(j+Ntsteps)]
m_vect.append(int(m[0]))
for i in range (0,2):
t = np.linspace(0, 1000, num=m[i]+1) # DO NOT CHANGE THIS PART
dt=(t[-1]-t[0])/(m[i]+1)
y1 = np.zeros(m[i]+1)
y2 = np.zeros(m[i]+1)
# Set the initial condition
y1[0] = y0[0]
y2[0]= y0[1]
for k in range(0,m[i]):
yk = np.array([y1[k],y2[k]])
yk1 = BEM(yk,f,L,M,dt)
y1[k+1]=yk1[0]
y2[k+1]=yk1[1]
T.append(t) #Save the time stepping
Y1.append(y1) #Save the Y1 values for m and 2m
Y2.append(y2) #Save the Y2 values for m and 2m
#################################
# RMSE PART 2 #
#################################
Y1_m=interpolate.interp1d(T[index],Y1[index], kind='linear')
Y1_m_Size2m=Y1_m(T[index+1])
Y2_m=interpolate.interp1d(T[index],Y2[index], kind='linear')
Y2_m_Size2m=Y2_m(T[index+1])
Y_m=np.array([[Y1_m_Size2m],[Y2_m_Size2m]])
Y_2m=np.array([[Y1[index+1]],[Y2[index+1]]])
rmse_Y_m=np.linalg.norm(Y_m-Y_2m)/(np.sqrt(2*m[0]))
RMSE_Y_m.append(float(rmse_Y_m))
index+=2
#################################
# PLOT PART 2 #
#################################
plt.figure("Root mean squared error rmse(m)\nPART 2")
plt.title("Root mean squared error")
plt.xlabel('Number of steps m $(5k, 7k, 9k, 11k, 13k, 15k)$')#Set a label to horizontal axis
plt.ylabel('rmse(m)')#Set a label to vertical axis
plt.grid()
plt.loglog(m_vect,RMSE_Y_m,color="orange",linewidth=2);
Code part 3
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#################################
# INITIALISATION #
#################################
#COEFFICIENT ###
## 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
y0 = np.array([1.060660171777250 +10, 2.121320343560527e+02])
#################################
# Solve with SciPy's function #
#################################
def system(t,state,a,g,b):
y1,y2=state
dy1_dt=a1*y1**g11*y2**g21-b1*y1
dy2_dt=a2*y1**g12*y2**g22-b2*y2
return [dy1_dt,dy2_dt]
p=([a1,a2],[[g11,g12],[g21,g22]],[b1,b2])
time=[0,1000]
res_ex=solve_ivp(system,time,y0,method="LSODA",args=p,atol=0.000001)
#################################
# PLOT PART 3 #
#################################
fig,ax=plt.subplots()
plt.title("Evolution in time of $y_1$ and $y_2$\nPART 3")
ax.plot(res_ex.t,res_ex.y[0,:],"-r")
ax.set_xlabel("Time $t$ [s]",fontsize=14)
ax.set_ylabel("$y_1$",color="red",fontsize=14)
ax.tick_params(axis="y",colors='red')
ax2=ax.twinx()
ax2.plot(res_ex.t,res_ex.y[1,:],"-b")
ax2.set_ylabel("$y_2$",color="blue",fontsize=14)
ax2.tick_params(axis="y",colors='blue')