explicit and implicit

AdamsWilliamson.py

+from cmath import pi
 import numpy as np
 import matplotlib.pyplot as plt
 import csv
@@ -23,20 +24,20 @@ def read_data(file):
     return np.column_stack([depth]), np.column_stack([rho]), np.column_stack([vp]), np.column_stack([vs])
 
 
-def mass(r, rho_r):
+def mass(rho, r):
     '''Calculates mass of earth inside sphere of radius r.'''
-    return (4/3)*np.pi*rho_r*r**3
+    return (4/3)*pi*rho*r**3
 
 
 # Berechnung der Dichte
-def density_expansion(rho, vp, vs, r1, r2):
+def density_expansion(rho_old, vp, vs, r_new, r_old):
     '''Returns Euler approximation for rho at r2.
     Method: 
     rho(r+h) = rho(r) + h*d(rho)/dr(r)'''
     G = 6.67 * 10**(-11)  # Gravitationskonstante
     phi = vp**2 - (4/3) * vs**2
-    rho_r = rho + (r2-r1) * (-G * mass(r1, rho) * rho / (r1**2*phi))
-    return rho_r
+    rho_new = rho_old + (r_new - r_old) * (-G * mass(r_old, rho_old) * rho_old / (r_old**2 * phi))
+    return rho_new
 
 
 file = os.getcwd() + r"\VpVs-valuesak135.csv"   
@@ -52,13 +53,23 @@ rho_0 = 3000  # kg/m^3
 M_0 = 5.973 * 10**24  # Masse in 18km Tiefe
 
 rho = np.zeros(len(density))
-rho[8] = rho_0
+rho[:9] = rho_0
 
 for i in range(9,len(rho)):
-    rho[i] = density_expansion(rho[i-1], vp[i-1], vs[i-1],r[i-1], r[i])
+    rho[i] = density_expansion(rho[i-1], vp[i-1], vs[i-1],r[i], r[i-1])
 
 
-print(rho)
+def corrected_mass(rho, r):
+    res = 0
+    for i in range(len(rho)):
+        res += (4/3)*pi*rho[i]*(r[i]**3 - r[i-1]**3)
+    return res
+
+print(corrected_mass(rho,r))
+
+plt.plot(rho, r)
+plt.plot(density, r)
+plt.show()
 
 # Berechnung von rho und M für jede Schicht
 # for i in range(8, len(rho) - 1):  # Die ersten Werte werden weg gelassen

AdamsWilliamsonCMB.py

+from cmath import pi
 import numpy as np
 import matplotlib.pyplot as plt
 import csv
@@ -25,13 +26,11 @@ def read_data(file):
 
 def mass(r, rho_r):
     '''Calculates mass of earth inside sphere of radius r.'''
-    mass = (4/3)*np.pi*rho_r*r**3
-    return mass
+    return (4/3)*np.pi*rho_r*r**3
 
 
-# Berechnung der Dichte
-def density_expansion(rho, vp, vs, r1, r2):
-    '''Returns Euler approximation for rho at r2.
+def density_expansion_explicit(rho, vp, vs, r1, r2):
+    '''Returns explicit Euler approximation for rho at r2.
     Method: 
     rho(r+h) = rho(r) + h*d(rho)/dr(r)'''
     G = 6.67 * 10**(-11)  # Gravitationskonstante
@@ -40,6 +39,23 @@ def density_expansion(rho, vp, vs, r1, r2):
     return rho_r
 
 
+def density_expansion_implicit(rho, vp, vs, r1, r2):
+    '''Returns implicit Euler approximation (Heun method) for rho at r2.
+    Method: 
+    rho(r+h) = rho(r) + h*d(rho)/dr(r+1) with d(rho)/dr(r+1) = d(rho)/dr(rho(r) + h*d(rho)/dr(r))'''
+    G = 6.67 * 10**(-11)  # Gravitationskonstante
+    phi = vp**2 - (4/3) * vs**2
+    rho_r = rho + (r2-r1) * (-G * mass(r2, density_expansion_explicit(rho, vp, vs, r1, r2)) * density_expansion_explicit(rho, vp, vs, r1, r2) / (r1**2*phi))
+    return rho_r
+
+def mass_correction(rho, r):
+    res = 0
+    for i in range(1,len(r)):
+        res += (4/3) * pi * rho[i] * -(r[i]**3-r[i-1]**3)
+    return res
+
+
+
 file = os.getcwd() + r"\VpVs-valuesak135.csv"   
 
 # Einlesen der Daten
@@ -52,20 +68,33 @@ r = 6371 - depth
 rho_0 = 3000  # kg/m^3
 M_0 = 5.973 * 10**24  # Masse in 18km Tiefe
 
-rho = np.zeros(len(density))
+rho_ex = np.zeros(len(density))
+rho_im = np.zeros(len(density))
 masse = np.zeros(len(density))
 masse[8] = M_0
-rho[8] = rho_0
+rho_ex[8] = rho_0
+rho_im[8] = rho_0
 
-for i in range(9, len(rho)):
+for i in range(9, len(rho_ex)):
     if depth[i] == 2891.5:
-        rho[i] = 9400
+        rho_ex[i] = 9900
+        rho_im[i] = 9900
     else:
-        rho[i] = density_expansion(rho[i-1], vp[i-1], vs[i-1],r[i-1], r[i])
-        masse[i] = mass(r[i-1], rho[i-1])
+        rho_ex[i] = density_expansion_explicit(rho_ex[i-1], vp[i-1], vs[i-1], r[i-1], r[i])
+        rho_im[i] = density_expansion_implicit(rho_im[i-1], vp[i-1], vs[i-1], r[i-1], r[i])
+        masse[i] = mass(r[i-1], rho_ex[i-1])
+
+
+print("Earth mass explicit: ", mass_correction(rho_ex, r))
+print("Earth mass implicit: ", mass_correction(rho_im, r))
+print("Earth mass exact: ", mass_correction(density, r))
+
 
-print(rho)
-print(masse)
+plt.plot(depth,rho_ex)
+plt.plot(depth,rho_im)
+plt.plot(depth,density)
+plt.legend(["ex", "im", "exact"])
+plt.show()
 
 # Berechnung von rho und M für jede Schicht
 # for i in range(8, len(rho) - 1):  # Die ersten Werte werden weg gelassen