Report -Bone remodeling

################################# # 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]

################################# # 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

################################# # 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)

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

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);

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')