text:
Tutorials

# Notebook - Mandyoc tutorial

Authors
Agustina Pesce
Victor Sacek
Jamison Assunção
Rafael Monteiro da Silva

# #Mandyoc: MANtle DYnamics simulatOr Code

Mandyoc is a finite element code written on top of the PETSc library to simulate thermo-chemical convection of the Earth’s mantle. Different linear and non-linear rheologies can be adopted, appropriately simulating the strain and stress pattern in the Earth’s crust and mantle, both in extensional or collisional tectonics. Additionally, the code allows variations of boundary condition for the velocity field in space and time, simulating different pulses of tectonism in the same numerical scenario

To simulate mantle thermochemical convection, we adopted the formulation for non-Newtonian fluids together with the Boussinesq approximation to solve the following equations of conservation of mass, momentum and energy:

$u_{i,i} = 0$
(1)#
$\sigma_{i j,j} + g_{i} \rho_{0}(1 - \alpha(T - T_{0})) = 0$
(2)#
$\frac{\partial T}{\partial t} + u_{i} T_{,i} = \kappa T_{,ii} + \frac{H}{c_p} - \frac{\alpha T g u_{e}}{c_{p}}$
(3)#

where $u$ is the velocity in the $i$ direction, $g$ is the gravity acceleration, $\rho_{0}$ is the reference rock density at temperature $T_0$, $\alpha$ is the coefficient of volumetric expansion, $T$ is the temperature, $\kappa$ is the thermal diffusivity, $H$ is the rate of radiogenic heat per unit of mass, $c_{p}$ is the specific heat, $\delta_{ij}$ is the Kronecker delta, and $\sigma_{ij}$ is the stress tensor:

$\sigma_{i j} = -P \delta_{i j} + \eta (u_{i,j} + u_{j,i})$
(4)#

where $P$ is the dynamic pressure and $\eta$ is the effective viscosity of the rock.

The equations of conservation of mass, momentum and energy are solved using the finite element method. Mandyoc uses quadrilateral elements for two dimensional grids. The mass and momentum equations subsection presents the numerical methods used to solve the mass and momentum equations and the energy equation subsection shows how an implicit formulation was used to solve the energy equation.

To simulate any scenario, the user must provide the parameter file param.txt and, if set, the necessary ASCII files with the initial temperature field, velocity field and/or the initial interfaces of the model. The flowchart showing the steps Mandyoc takes to solve the equations of conservation of mass, momentum and energy.

It shows that once the code starts running and the input files are read, Mandyoc uses the effective viscosity field $\eta$ to calculate the velocity field $u$ and checks if the convergence condition satisfies the tolerance $tol$.

While the minimum tolerance is not reached, Mandyoc utilizes the Uzawa's method to iteratively calculate a new $u$ and $P$ fields. The updated fields modify the viscosity field $\eta$, which in turn disturbs the velocity field again. These fields are updated until tolerance is reached.

Additionally, the compositional factor $C$ is evaluated for an advection as:

$\frac{\partial C}{\partial t} + u_{i}C_{,i} = 0$
(5)#

Its solution is calculated placing randomly a number of particles within each finite element of the mesh, which are displaced based on the adjacent node velocity values. The individual value for each particle is obtained by linear interpolation of the node values.

List of scientific works done with Mandyoc or its earlier version:

# #Showcase: Continental rift

This example simulates the evolution of divergent margins, taking into account the plastic rheology and the sin-rift geodynamics. The idea of this tutorial is to explain how to generate the initial model to reach the following result:

## #Input files

There are a few ways to run this model. The one that will present comes from the assumption that the geometry, temperature and velocity at the beginning of the simulation are know in a reasonable fashion. Therefore, we must provide them as input files, which are:

• an interface file with the boundary of the different layers and their composition information,
• an initial temperature file with the initial temperature field,
• and an initial velocity file with the initial velocity field

These files describe the following setup model for the continental rift:

# Python packages:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

# Import wornings package to ignore any wrnings:
import warnings
warnings.filterwarnings('ignore')

# Import some functions that I created for this tutorial:
from functions import (
diffusion_equation,
colors
)

The domain of the model comprises 1600 x 300 km2 that is discretized by elements of 1 x 1 km2:

lx, lz = 1600.0e3, 300.0e3

nx, nz = 1601, 301

The regular meshgrid can be easily created using:

x = np.linspace(0, lx, nx)
z = np.linspace(-lz, 0, nz)

xx, zz = np.meshgrid(x, z)

### #Interface file

For this simulation, the domain is subdivided into a few lithologial units:

• asthenosphere,
• lithospheric mantle,
• weak seed,
• lower crust,
• upper crust, and
• sticky air.

Each of them with their different compositional properties.

Additionally, to ensure the nucleation of rifting at the center of the numerical domain, a weak seed will add in the lithospheric mantle.

h_sair = 40.0e3

h_lower_crust = 20.0e3

h_upper_crust = 20.0e3

h_litho = 130.0e3

# Seed depth position under the crust:
z_seed = 13.0e3

The different lithologies can be expressed in a dictionary:

interfaces = {
"litho": -np.ones(nx) * h_litho,
"seed_base": -np.ones(nx) * (z_seed + h_lower_crust + h_upper_crust),
"seed_top": -np.ones(nx) * (z_seed + h_lower_crust + h_upper_crust),
"lower_crust": -np.ones(nx) * (h_lower_crust + h_upper_crust),
"upper_crust": -np.ones(nx) * h_upper_crust,
"air": np.zeros(nx),
}
for i in interfaces:
interfaces[i] = interfaces[i] - h_sair

The considerations above result in the geometry shown in the following plot:

fig, ax = plt.subplots(figsize=(16, 8))

for label, layer in interfaces.items():
ax.plot(x/1e3, (layer + h_sair)/1e3, label=f"{label}")

ax.set_yticks(np.arange((-lz + h_sair)/1e3, h_sair/1e3, 10))
ax.set_xlim([0, lx/1e3])
ax.set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
ax.set_xlabel("x [km]")
ax.set_ylabel("Depth [km]")

plt.title("Interfaces")
plt.legend(loc='upper right')

plt.show()

The weak seed can be created:

# Seed horizontal position:
x_seed = 750.0e3
# Seed thickness:
h_seed = 6.0e3
d_x = lx / (nx - 1)

cond = (x_seed - d_x <= x) & (x < x_seed + d_x)

# Add a thickness in the seed layer:
interfaces["seed_base"][cond] = interfaces["seed_base"][cond] - h_seed//2

interfaces["seed_top"][cond] = interfaces["seed_top"][cond] + h_seed//2

Plot again the interfaces:

fig, ax = plt.subplots(figsize=(16, 8))

for label, layer in interfaces.items():
ax.plot(x/1e3, (layer + h_sair)/1e3, label=f"{label}")

ax.set_yticks(np.arange((-lz + h_sair)/1e3, h_sair/1e3, 10))
ax.set_xlim([0, lx/1e3])
ax.set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
ax.set_xlabel("x [km]")
ax.set_ylabel("Depth [km]")

plt.title("Interfaces")
plt.legend(loc='upper right')

plt.show()

Now we can create the interface file. This file contains the layers properties and the interface's depth between these layers.

We are going to adopt a visco-plastic rheology where the effective viscosity $\eta$ combines nonlinear power law viscous rheology and a plastic yield criterion. The viscous component is given by:

$\eta_{visc} = C A^{-1/n} \dot{\varepsilon}^{(1-n)/n} \exp{\frac{Q + V P}{n R T}}$
(6)#

So, the layers properties that we need to specify are:

• Compositional factor ($C$)
• Density ($rho$)
• Radiogenic heat ($H$)
• Pre-exponential scale factor ($A$)
• Power law exponent ($n$)
• Activation energy ($Q$)
• Activation volume ($V$)

To simulate the free surface, the sticky air approach is adopted, taking into account a 40 km thick layer with a relatively low viscosity material but with a compatible density with the atmospheric air.

AsthenosphereLitho mantleSeedLitho mantle over the seedLower crustUpper crustAir
C-1.01.00.11.01.01.01.0
rho
$kg/m^3$
(7)#
3378.03354.03354.03354.02800.02700.01.0
H
$W/kg$
(8)#
0.09.0e-129.0e-129.0e-12HlcHuc0.0
A
$Pa^{-n}/s$
(9)#
1.393e-142.4168e-152.4168e-152.4168e-158.574e-288.574e-281.0e-18
n-3.03.53.53.54.04.01.0
Q
$kJ/mol$
(10)#
429.0e3540.0e3540.0e3540.0e3222.0e3222.0e30.0
V
$m^3/mol$
(11)#
15.0e-625.0e-625.0e-625.0e-60.00.00.0
# Define the density and the radiogenic heat for the upper and lower crust:
rho_uc = 2700.0
rho_lc = 2800.0
Huc = 2.5e-6 / rho_uc
Hlc = 0.8e-6 / rho_lc
# Define the layer properties:
layer_properties = f"""
C   1.0       1.0        0.1        1.0        1.0         1.0         1.0
rho 3378.0    3354.0     3354.0     3354.0     {rho_lc}    {rho_uc}    1.0
H   0.0       9.0e-12    9.0e-12    9.0e-12    {Hlc}       {Huc}       0.0
A   1.393e-14 2.4168e-15 2.4168e-15 2.4168e-15 8.574e-28   8.574e-28   1.0e-18
n   3.0       3.5        3.5        3.5        4.0         4.0         1.0
Q   429.0e3   540.0e3    540.0e3    540.0e3    222.0e3     222.0e3     0.0
V   15.0e-6   25.0e-6    25.0e-6    25.0e-6    0.0         0.0         0.0
"""
# Create and save the interface file:
with open("interfaces.txt", "w") as f:
# Write the layer properties
for line in layer_properties.split("\n"):
line = line.strip()
if len(line):
f.write(" ".join(line.split()) + "\n")

# Write the interfaces of the layers adding the thickness of the sticky air layer:
data = np.array(tuple(interfaces.values())).T
np.savetxt(f, data, fmt="%.1f")

### #Initial temperature file

The initial temperature structure is depth dependent and is 0°C at the surface and 1300°C at the base of the lithosphere at 130 km. With these boundary conditions, the initial temperature structure in the interior of the lithosphere is given by the solution of 1-D diffusion equation:

$\kappa \frac{\partial^2 T(z)}{\partial z^2} + \frac{H(z)}{c_p} = 0$
(12)#

where $H(z)$ is the internal heat production of the different layers.

The sublithospheric temperature follows an adiabatic increase up to the bottom of the model:

$T = T_p \exp{\frac{g \alpha z}{c_p}}$
(13)#

Where $T_p$ is the potential temperature for the mantle, $g$ is the gravity value, $\alpha$ is the volumetric expansion coefficient, $c_p$ is the specific heat capacity.

• $\kappa = 1x10^{-6}$ $m²/s$
• $c_p = 1250$ $J kg^{-1} K^{-1}$
• $T_p = 1262$ °C
• $g = -10$ $m/s²$
• $\alpha = 3.28x10^{-5}$ $K^{-1}$
kappa = 1.0e-6
ccapacity = 1250
tem_p = 1262
g = -10
alpha = 3.28e-5
# Temperature when z is over 130 km:
temp_z = 1300 * (-z - h_sair) / (h_litho)

# Sublithospheric temperature:
temp_adiabatic = tem_p / np.exp(g * alpha * (-z - h_sair) / ccapacity)

# T = 0 when z is over zero:
temp_z[z + h_sair >= 0.0] = 0.0

# When z is below 130 km (Lithosperic depth), the T = temp_adiabatic:
temp_z[temp_z > temp_adiabatic] = temp_adiabatic[temp_z > temp_adiabatic]

Now, we will apply the 1-D diffusion equation in the model.

First, create the internal heat production model:

# Initialize H as a grid of zeros:
H = np.zeros_like(temp_z)

# Add the H value for the upper crust:
cond = (z <= interfaces["air"][0]) & (z > interfaces["upper_crust"][0])
H[cond] = Huc

# Add the H value for the lower crust:
cond = (z <= interfaces["upper_crust"][0]) & (z > interfaces["lower_crust"][0])
H[cond] = Hlc

Apply the diffusion equation for 10000 years:

dt = 10000
dz = lz / (nz - 1)

# Apply the 1-D diffusion equation:
temp_z = diffusion_equation(
temp_z,
z,
interfaces,
H,
dt,
dz,
kappa,
ccapacity
)

Save the temperature field in a file:

temp_z = np.ones_like(xx) * temp_z[:, None]

print(np.shape(temp_z))

# Save the initial temperature file
np.savetxt(
"input_temperature_0.txt",
np.reshape(temp_z, (nx * nz)),
)

Plot the temperature field and a thermal profile:

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 8), sharey=True)

# Plot temperature field:
img_temp = axs[0].contourf(
xx/1.0e3,
(zz + h_sair)/1.0e3,
temp_z,
levels=np.arange(0, 1610, 100)
)
cbar = fig.colorbar(img_temp, orientation='vertical', ax=axs[0])
cbar.set_label("Temperature [°C]")

axs[0].set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
axs[0].set_ylabel("Depth [km]")
axs[0].set_xlabel("x [km]")

# Plot temperature profile:
axs[1].set_title("Thermal profile [°C]")
axs[1].plot(temp_z[:, 0], (z + h_sair)/1.0e3, "-k")
for label, layer in interfaces.items():
axs[1].plot(x/1e3, (layer + h_sair)/1e3, label=f"{label}")

axs[1].set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
axs[1].set_xlabel("Temperature [°C]")
axs[1].set_ylabel("Depth [km]")
plt.legend(loc="lower left")

plt.show()

### #Initial velocity field

The velocity field simulate the lithospheric stretching assuming a reference frame fixed on the lithospheric plate on the left side of the model, and the plate on the right side moves rightward with a velocity of $v_{rift} = 1$ cm/year. The velocity field in the left and right boundaries of the model is chosen to ensure conservation of mass and is symmetrical if the adopted reference frame movies to the right with a velocity of $v_{rift}/2$ cm/year relative to the left plate. The inflow on the vertical boundaries below the lithosphere compensates for outflow.

So, the horizontal velocity field along the left and right borders of the domain are:

$velocity_x(x = 0) = \left\{ \begin{array}{r l} 0, & \text{z > -h_c - h_{sair}} \\ -v \frac{z + h_c + h_{sair}}{h_a}, & \text{z \leq -h_c - h_{sair}} \end{array} \right.$
(14)#
$velocity_x(x = lx) = \left\{ \begin{array}{rl} 0, & \text{z > h_{sair}} \\ v_{rift}, & \text{z > -h_c - h_{sair}} \\ v \frac{z + h_c + h_{sair}}{h_a} + v_{rift}, & \text{z \leq -h_c - h_{sair}} \end{array} \right.$
(15)#

where $v = v_{rift} \frac{lz - h_{sair}}{h_a}$ so that:

$\int^{0}_{-lz} velocity_x(x = 0) = \int^{0}_{-lz} velocity_x(x = lx) = \frac{v h_a}{2}$
(16)#

$h_c = 150$ km is the thickness of the upper layer with constant velocity, corresponding to the lithosphere $h_{litho} = 130$ km and 20 km of the asthenosphere and $h_a = lz - h_c - h_{sair}$

# Converting 1 cm/year to m/s:
v_rift = 1e-2 / (365 * 24 * 3600)
h_c= 150.0e3
h_a = lz - h_sair - h_c
v = v_rift * (lz - h_sair) / h_a
# Creating horizontal and vertical velocity:
velocity_xx, velocity_zz = np.zeros_like(xx), np.zeros_like(xx)

Velocity for the left side (x == 0):

cond = (xx == 0) & (zz <= -h_c - h_sair)
velocity_xx[cond] += -v * (zz[cond] + h_c + h_sair) / h_a

Velocity for the rigth side (x == lx):

cond = (xx == lx) & (zz > -h_c - h_sair)
velocity_xx[cond] += v_rift

cond = (xx == lx)& (zz <= -h_c - h_sair)
velocity_xx[cond] += v * (zz[cond] + h_c + h_sair) / h_a + v_rift

# Consider 10 km of the air layer with constant velocity:
frac_air = 10e3
velocity_xx[zz >= -h_sair + frac_air] = 0

Due to the mass conservation is assumed, the sum of the integrals over the boundaries (material flow) must be zero.

# For the left side:
v0 = velocity_xx[xx == 0]
sum_velocity_left = np.sum(v0[1:-1]) + (v0[0] + v0[-1]) / 2.0

# For the right side:
vf = velocity_xx[(xx == lx)]
sum_velocity_right = np.sum(vf[1:-1]) + (vf[0] + vf[-1]) / 2.0

diff = (sum_velocity_right - sum_velocity_left) * dz

print("Sum of the integrals over the boundary is:", diff)

When the sum of the integrals on the boundaries is not zero due to rounding errors, we add the difference back on the top to compensate make it zero. Despite not being negligible, this is reasonable because it is a very small correction.

velocity_zz[zz == 0] = -diff / lx

Plot the velocity:

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(9, 9), sharey=True)

# Plot horizontal component:
axs[0].plot(velocity_xx[:, 0], (z + h_sair)/1e3, "k-", label="left side")
axs[0].plot(velocity_xx[:, -1], (z + h_sair)/1e3, "r-", label="right side")

axs[0].legend()
axs[0].grid()
axs[0].set_ylim([(-lz + h_sair)/1e3 , h_sair/1e3])
axs[0].set_xlim([-8e-10, 8e-10])
axs[0].set_xlabel(" Velocity [m/s]")
axs[0].set_ylabel("Depth [km]")
axs[0].set_title("Horizontal component of velocity")

# Plot vertical component
axs[1].plot(velocity_zz[:, 0], (z + h_sair)/1e3, "k-", label="left side")
axs[1].plot(velocity_zz[:, -1], (z + h_sair)/1e3, "r-", label="right side")

axs[1].legend()
axs[1].grid()
axs[1].set_xlim([-8e-10, 8e-10])
axs[1].set_xlabel(" Velocity [m/s]")
axs[1].set_title("Vertical component of velocity")

plt.show()

Create and save the initial velocity file:

velocity = np.zeros((2, nx * nz))

velocity[0, :] = np.copy(np.reshape(velocity_xx, nx * nz))
velocity[1, :] = np.copy(np.reshape(velocity_zz, nx * nz))

velocity = np.reshape(velocity.T, (np.size(velocity)))

np.savetxt("input_velocity_0.txt", velocity, header="v1\nv2\nv3\nv4")

### #Create the parameter file

The parameter file contains the information that is necessary for the simulation to run. There are many parameters to define, so we will only show/discuss some of the important ones for this simulation. But I invite you to visit the Mandyoc documentation website to know more about the parameter file.

• Define the number of Lagrangian particles in each element, particles_per_element = 40.

• Specify the maximum time and step of the simulation with step_max = 5000 and time_max = 100.0e6 (100 My).

• For the temperature boundary conditions, the temperature values is fixed on all sides during the simulation, with the nodes of the model in the air maintained at 0°C.

• For the velocity boundary conditions, free slip is applied on the top and bottom of the numerical domain.

• To avoid a symmetric, laterally homogeneous model, we introduced a random perturbation of the initial strain in each finite element of the model.

• Set the pre-defined rheology model to use during simulation. In this case: rheology_model = 9

params = f"""
# Geometry
nx                                  = {nx}          # n. of nodes in the horizontal direction
nz                                  = {nz}          # n. of nodes in the vertical direction
lx                                  = {lx}          # extent in the horizontal direction
lz                                  = {lz}          # extent in the vertical direction

# Simulation options
denok                               = 1.0e-15       # tolerance criterion for the Uzawa's scheme, default is 1.0e-4
Xi_min                              = 1.0e-7        # tolerance criterion for the convergence of the non-linear flow, default is 1.0e-14
random_initial_strain               = 0.3           # non-dimensional value for the initial strain perturbation for the entire domain, default is 0.0
initial_dynamic_range               = True          # method to smoothen convergence of the velocity field in scenarios with wide viscosity range, default is False [True/False], Gerya (2019)

# Particles options
particles_per_element               = 40            # n. of Lagrangian particles in each element, default is 81
particles_perturb_factor            = 0.7           # indicates the amount of perturbation of the initial location of the particles relative to a regular grid distribution. Default is 0.5 [values are between 0 and 1]

# Time constrains
step_max                            = 5000          # maximum time-step of the simulation [steps]
time_max                            = 100.0e6       # maximum time of the simulation [years]
dt_max                              = 10.0e6        # maximum time between steps of the simulation [years]
step_print                          = 10            # make output files every <step_print>
sub_division_time_step              = 0.5           # re-scale value for the calculated time-step, default is 1.0

# Viscosity
viscosity_reference                 = 1.0e26        # reference mantle viscosity [Pa.s]
viscosity_max                       = 1.0e25        # maximum viscosity during simulation [Pa.s]
viscosity_min                       = 1.0e18        # minimum viscosity during simulation [Pa.s]
viscosity_per_element               = constant      # sets if viscosity is constant or linearly variable for every element, default is variable [constant/variable]
viscosity_mean_method               = arithmetic    # defines method do calculate the viscosity for each element, default is harmonic [harmonic/arithmetic]
viscosity_dependence                = pressure      # defines if viscosity depends on pressure or depth, default is depth [pressure/depth]

# External ASCII inputs/outputs
interfaces_from_ascii               = True          # set if interfaces between lithologies are read from an ASCII file (interfaces.txt), default is False [True/False]
n_interfaces                        = {len(interfaces.keys())}            # set the number of interfaces to be read from the interfaces ASCII file (interfaces.txt)
temperature_from_ascii              = True          # set if initial temperature is read from an ASCII file (input_temperature_0.txt), default is False [True/False]
velocity_from_ascii                 = True          # set if initial velocity field is read from an ASCII file (input_velocity_0.txt), default is False [True/False]
variable_bcv                        = False         # allows velocity field re-scaling through time according to an ASCII file (scale_bcv.txt), default is False [True/False]
multi_velocity                      = False         # set if boundary velocities can change with time from ASCII file(s) (multi_veloc.txt and additional input_velocity_[X].txt files), default is False [True/False]
binary_output                       = False         # set if output is in binary format, default is False [True/False]
print_step_files                    = True          # set if the particles position are printed to an output file, default is True [True/False]

sticky_blanket_air                  = True          # default is False [True/False]
precipitation_profile_from_ascii    = False         # default is False [True/False]
climate_change_from_ascii           = False         # default is False [True/False]

# Physical parameters
temperature_difference              = 1500.         # temperature difference between the top and bottom of the model (relevant if <temperature_from_ascii> is False) [K]
thermal_expansion_coefficient       = 3.28e-5       # value for the coefficient of thermal expansion [1/K]
thermal_diffusivity_coefficient     = 1.0e-6        # value for the coefficient of thermal diffusivity [m^2/s]
gravity_acceleration                = 10.0          # value for the gravity acceleration [m/s^2]
density_mantle                      = 3300.         # value for the mantle reference density [kg/m^3]
external_heat                       = 0.0e-12       # ok
heat_capacity                       = 1250.         # value for the heat capacity [J/K]

non_linear_method                   = on            # set if non linear method is used for the momentum equation [on/off]
adiabatic_component                 = on            # set if adiabatic heating/cooling is active [on/off]
radiogenic_component                = on            # set if radiogenic heating is active [on/off]

# Velocity boundary conditions
top_normal_velocity                 = fixed         # set the normal velocity on the top side of the model to be fixed or free [fixed/free]
top_tangential_velocity             = free          # set the tangential velocity on the top side of the model to be fixed or free [fixed/free]
bot_normal_velocity                 = fixed         # set the normal velocity on the bottom side of the model to be fixed or free [fixed/free]
bot_tangential_velocity             = free          # set the tangential velocity on the bot side of the model to be fixed or free [fixed/free]
left_normal_velocity                = fixed         # set the normal velocity on the left side of the model to be fixed or free [fixed/free]
left_tangential_velocity            = fixed         # set the tangential velocity on the left side of the model to be fixed or free [fixed/free]
right_normal_velocity               = fixed         # set the normal velocity on the right side of the model to be fixed or free [fixed/free]
right_tangential_velocity           = fixed         # set the tangential velocity on the right side of the model to be fixed or free [fixed/free]

surface_velocity                    = 0.0e-2        # ok

# Temperature boundary conditions
top_temperature                     = fixed         # set temperature on the top side of the model to be fixed or free [fixed/free]
bot_temperature                     = fixed         # set temperature on the bottom side of the model to be fixed or free [fixed/free]
left_temperature                    = fixed         # set temperature on the left side of the model to be fixed or free [fixed/free]
right_temperature                   = fixed         # set temperature on the right side of the model to be fixed or free [fixed/free]

rheology_model                      = 9             # flag n. of a pre-defined rheology model to use during simulation
"""

# Create the parameter file
with open("param.txt", "w") as f:
for line in params.split("\n"):
line = line.strip()
if len(line):
f.write(" ".join(line.split()) + "\n")

## #Run the model

In this example, mandyoc use the following flags:

-seed 0,2
-strain_seed 0.0,1.0

the first flag -seed specify which layers will have a specific initial strain. These layers are the number 0 and 2 for the asthenosphere and the mantle weak seed in the middle of the lithospheric mantle.

The second flag -strain_seed specify the magnitude of the initial strain, with respect to each layer indicated in the previous flag. In these case, the asthenosphere will have no initial strain and the mantle seed will have a initial constant strain of 1.0. The other layers will present a random perturbation for the value of the initial strain, and the maximum magnitude of this initial strain is defined by the parameter random_initial_strain in parameters file.

You can run the model as:

mpirun -n NUMBER_OF_CORES mandyoc -seed 0,2 -strain_seed 0.0,1.0

Where the NUMBER_OF_CORES value must be set by the user.

This process took around 2 days to finish the temporal evolution.

## #Post-processing

After Mandyoc finished running the model, different files are obtained with the temporal evolution of the model as it was specified in the parameter file. In this case, we are only going to use the density, cumulative strain, temperature, velocity, viscosity and density files.

### #Plot the results

Now, we will plot the Mandyoc output file for only 5 different steps times. These files are into output/ directory.

Determine the initial and final step to make the plots:

step_initial = 0
step_final = 5000

d_step = 1000

The dark and light orange colors represent the upper and lower crust, respectively, while the dark and light green colors represent lithospheric and sublithospheric mantle, respectively. Shades of gray indicate the magnitude of the cumulative strain. The blue and orange bars indicate the width of the extended continental crust in both conjugate margins.

for cont in range(step_initial, step_final + d_step, d_step):

# Read the files:
# Time
filename = "output/time_" + str(cont) + ".txt"
# Density
filename = "output/density_" + str(cont) + ".txt"
rho = load_data(filename, nx, nz)
# Strain
filename = "output/strain_" + str(cont) + ".txt"
strain = load_data(filename, nx, nz)
strain[rho < 200] = 0
strain_log = np.log10(strain)
# Temperature
filename = "output/temperature_" + str(cont) + ".txt"
temp = load_data(filename, nx, nz)
# Velocity
filename = "output/velocity_" + str(cont) + ".txt"
veloc_x, veloc_z = load_velocity(filename, nx, nz)
# Density
filename = "output/viscosity_" + str(cont) + ".txt"
visco = load_data(filename, nx, nz)
visco_log = np.log10(visco)
visco_log[visco_log == -np.inf] = 0

# Create the plots:
fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(20, 10), sharex=True)
axs[0].set_title("Time = %.1lf Myr\n\n" % (time[0]/1.0e6))

# Plot density and strain_log
# Density
img_rho = axs[0].contourf(
xx/1e3,
(zz + h_sair)/1e3,
rho,
levels=[100, 2700, 2800, 3354, 3378],
colors=colors,
)
plt.colorbar(
img_rho,
ax=axs[0],
orientation="vertical",
label="Density [$kg/m^3$]"
)
# Strain_log
strain_aux = strain_log.copy()
strain_aux[strain_aux == -np.Inf] = np.nan
strain_aux[strain_aux < -0.5] = np.nan
strain_aux[strain_aux > 0.9] = np.nan
img_strain = axs[0].scatter(
xx/1e3,
(zz + h_sair)/1e3,
c=strain_aux,
zorder=100,
s=0.5,
alpha=0.05,
cmap=plt.get_cmap("Greys"),
vmin=-0.5,
vmax=0.9,
)
axs[0].axis('equal')
axs[0].set_ylabel("Depth [km]")

# Plot temperature and velocity
# Temperature
img_temp = axs[1].pcolormesh(
xx/1e3,
(zz+h_sair)/1e3,
temp,
cmap="viridis",
)
plt.colorbar(
img_temp,
ax=axs[1],
orientation="vertical",
label="Temperature [°C]"
)
# Velocity
axs[1].quiver(
xx[::20, ::20]/1e3,
(zz[::20, ::20] + h_sair)/1e3,
(veloc_x[::20, ::20]/1e3).T,
(veloc_z[::20, ::20]).T/1e3,
color="w"
)
axs[1].axis('equal')
axs[1].set_ylabel("Depth [km]")

# Plot viscosity
img_visco = axs[2].pcolormesh(
xx/1e3,
(zz+h_sair)/1e3,
visco,
cmap="viridis",
norm=LogNorm(),
)
plt.colorbar(
img_visco,
ax=axs[2],
orientation="vertical",
label="Viscosity [Pa s]"
)
axs[2].axis('equal')
axs[2].set_ylabel("Depth [km]")
axs[2].set_xlabel("x [km]")

plt.show()


text:
# An orientation on a $(n-1)$-dimensional submanifold.

This question goes on where this question ended. Given is an non-empty $(n-1)$-dimensional ($n\ge2$) differentiable submanifold $X\subset\mathbb{R}^n$ such that there exists an open $U\subset\mathbb{R}^n$ with $X\subset U$ and a differentiable function $f:U\to\mathbb{R}$ such that $X=f^{-1}(0)$ and $Df(p)\neq0$ for all $p\in X$. Then $X$ is (as shown in the other question) the topological boundary in $U$ of $A:=\{p\in U\,|\,f(p)<0\}$.

I have to show that $X$ is an orientable manifold. That is, there is an differentiable atlas $\mathcal{A}$ such that $\det(D(h_i\circ h_j^{-1})(h_j(p))>0$ for all $p\in U_i\cap U_j$, for all charts $(U_i,h_,V_i),(U_j,h_j,V_j)\in\mathcal{A}$. Alternatively, there exists a $(n-1)$-form on $X$ which vanishes nowhere.

Now intuitively I feel like this is clear. Let $V\subset X$ be a connected component; I think that it is the case that $X$ is actually the smooth boundary of $A$, since if it wasn't, then $V$ (or some other connected component) would have some sort of kink somewhere, which would contradict it being a submanifold; as $A$ is open in $\mathbb{R}^n$ (since $U$ is open and $A$ is open in $U$), $A$ is orientable (restricting the standard orientation on $\mathbb{R}^n$ to $A$), and an orientation on $A$ induces an orientation on $\partial A=X$. Can this be turned into a formal argument?

Alternatively, how would one approach a proof using the two definitions of orientation given above? As the gradient $Df(p)$ is non-zero everywhere on $X$, and $X$ is the level set of $f$ at $C=0$, the gradient would be an ideal candidate to provide the outward pointing normal vector needed.

Any help is much appreciated. This question has kept me busy for some time now.

• It would be helpful if you could explain what $A$ is without referring to some answer to some other question. – Thomas May 10 '16 at 19:03
• You are right, I have edited. It wasn't in the answer of some other question, but in the question itself, though. – B. Pasternak May 10 '16 at 19:04

The following $n-1$ form $\omega$ on $X$ is nowhere zero and thus orients $X$:
For $x\in X$ and $v_1,v_2,\cdots,v_{n-1}\in T_x(X)=(\operatorname {grad}f(x) )^\perp\subset T_x(\mathbb R^n)= \mathbb R^n$ define $$\omega_x(v_1,v_2,\cdots,v_{n-1})=\det (v_1,\cdots,v_{n-1},\operatorname {grad}f(x) ).$$

format: plain


text:
### Home > AC > Chapter 5 > Lesson 5.1.6 > Problem5-59

5-59.
1. Consider the rule y = x − 2 . 5-59 HW eTool (Desmos). Homework Help ✎

1. Without graphing, find the zero of y = x − 2. Remember that a zero of a rule is an input that makes the output zero.

2. Make a table and graph y = x − 2 on graph paper.

3. How could you find the zero of y = x − 2 with your graph from part (b)? What about with the table? Explain.

Set the equation equal to zero, since your output would be zero, and then solve for x.

Find the point where y = 0.

Use the eTool below to help you with this problem.
Click on the link at right for the full eTool version: AC 5-59 HW eTool


text:
Board index Contests Integration Contest - Season 4

Integration Contest - Season 4

This is not an ordinary forum. This is the place for the most brilliant ideas and techniques.

Moderators: Random Variable, sos440

Re: Integration Contest - Season 4

Tue Mar 11, 2014 9:03 pm
galactus
Global Moderator

Posts: 902
Well, RV, I decided to try this one using this method. I will probably screw it up, but see if it's OK.

$\displaystyle \int_{0}^{1}\frac{\sqrt{x(x-1)}}{(x^{2}+1)^{2}}dx$

The poles are at $\displaystyle \pm i$

order 2, so diffing is in order.

$\displaystyle 2\pi i Res(i)=\frac{-1}{8\sqrt{i-1}}$

$\displaystyle 2\pi i Res(-i)=\frac{-1}{8\sqrt{-1-i}}$

add them gives $\displaystyle \frac{-\pi}{4}\sqrt{\sqrt{2}-1}i$

the residue at infinity is 0.

I tried the same args.

below the cut: $\displaystyle arg(z)=2\pi, \;\ arg(z-1)=\pi$

above the cut: $\displaystyle arg(z)=0, \;\ arg(z-1)=\pi$

$\displaystyle \int_{0}^{1}\frac{x^{1/2}e^{2\pi i/2}(x-1)^{1/2}e^{\pi i/2}}{(x^{2}+1)^{2}}dx-\int_{0}^{1}\frac{x^{1/2}e^{i(0)}(x-1)^{1/2}e^{i\pi/2}}{(x^{2}+1)^{2}}dx=\frac{-\pi}{4}\sqrt{\sqrt{2}-1}i$

$\displaystyle -i\int_{0}^{1}\frac{x^{1/2}(1-x)^{1/2}}{(x^{2}+1)^{2}}dx-i\int_{0}^{1}\frac{x^{1/2}(x-1)^{1/2}}{(x^{2}+1)^{2}}dx=\frac{-\pi}{4}\sqrt{\sqrt{2}-1}i$

$\displaystyle -2i\int_{0}^{1}\frac{x^{1/2}(x-1)^{1/2}}{(x^{2}+1)^{2}}dx=\frac{-\pi}{4}\sqrt{\sqrt{2}-1}i$

$\displaystyle \int_{0}^{1}\frac{\sqrt{x(x-1)}}{(x^{2}+1)^{2}}dx=\frac{\pi}{8}\sqrt{\sqrt{2}-1}$

Re: Integration Contest - Season 4

Thu Mar 13, 2014 1:37 pm
galactus
Global Moderator

Posts: 902
Here's another integral in case anyone is interested. Have we done this one before?. It sure seems familiar. I saw it on AoPS.

This one sure looks familiar, but I can find no reference to it. Strange.

$\displaystyle \int_{0}^{\frac{\pi}{2}}\frac{\sin^{2}(ax)}{\sin^{2}(x)}dx=\frac{\pi a}{2}+\frac{\sin(\pi a)}{2}+\frac{a\cdot \sin(\pi a)}{2}\left(\psi(a/2)-\psi\left(\frac{a+1}{2}\right)\right)$

when 'a' in an integer, then this has a general form of just $\displaystyle \frac{\pi a}{2}$. It is the non-integer values of 'a' that require more work and the formula involving digamma.

EDIT:
I made progress on this one by using the rather obscure identity:

$\displaystyle \int_{0}^{\pi}\sin^{p}(x)e^{iqx}dx=\frac{\pi}{2^{p}}\left[\frac{\Gamma(1+x)e^{\pi i q/2}}{\Gamma(1+\frac{p+q}{2})\Gamma(1+\frac{p-q}{2})}\right]...[1]$

I have to assemble the pieces though, but it's there. One may write:

$\displaystyle \int_{0}^{\frac{\pi}{2}}\frac{\sin^{2}(ax)}{\sin^{2}(x)}dx=1/2\int_{0}^{\pi}\frac{1-\cos(ax)}{1-\cos(x)}dx$

But first, we can evaluate $\displaystyle \int_{0}^{\pi}\frac{\sin(ax)}{\tan(x/2)}dx=\int_{0}^{\pi}\frac{\sin(ax)(1+\cos(x))}{\sin(x)}dx=\int_{0}^{\pi}\frac{\sin(ax)\sin(x)}{1-\cos(x)}dx=\int_{0}^{\pi}\frac{2\sin(ax)+\sin(a+1)x+\sin(a-1)x}{2\sin(x)}dx$

to be continued...............I have got the solution assembled. Please ask if you would like to see it. I had trouble with a nasty/tedious limit, but it appears to work out rather nicely. The integral identity in [1] can be shown by using contours and considering:

$\displaystyle f(x)=\frac{1}{(u+ix)^{a}(v-ix)^{b}}$ and make the sub $\displaystyle x=cot(t)$. This is one way to go about it.

EDIT:

I was playing around with this when u=v=1.

$\displaystyle \int_{-\infty}^{\infty}\frac{1}{(1+ix)^{a}(1-ix)^{b}}dx$

Factor out i:

$\displaystyle \frac{1}{i^{a}}\int_{-\infty}^{\infty}\frac{1}{(x-i)^{a}(1-ix)^{b}}dx$

There are poles at $\displaystyle x=\pm i$. Only $\displaystyle i$ is in the contour.

But, it has order 'a', so diff a-1 times to find the residue.

$\displaystyle =\frac{2\pi i}{i^{a}}\lim_{z\to i}\frac{1}{(a-1)!}\frac{d^{a-1}}{dz^{a-1}}\frac{1}{(1-iz)^{b}}$

do this a-1 times and note that $\displaystyle (a-1)!=\Gamma(a)$

$\displaystyle =2\pi\frac{b(b+1)(b+2)\cdot\cdot\cdot (b+a-2)}{\Gamma(a)}\lim_{z\to i}(1-iz)^{1-a-b}$

$\displaystyle =\pi 2^{2-a-b}\frac{b(b+1)\cdot\cdot\cdot (b+a-2)}{\Gamma(a)}$

$\displaystyle =2^{2-a-b}\pi \frac{\Gamma(a+b-1)}{\Gamma(a)\Gamma(b)}$

Re: Integration Contest - Season 4

Sat Mar 22, 2014 10:15 am
galactus
Global Moderator

Posts: 902
Evaluate

$\displaystyle \int_{0}^{\infty}\frac{\cos(ax)}{b^{2n}+x^{2n}}dx$

Re: Integration Contest - Season 4

Sat Mar 22, 2014 3:51 pm
Random Variable Integration Guru

Posts: 381
First we need to find the partial fraction decomposition of $\displaystyle \frac{1}{x^{2n}+b^{2n}}$.

$\displaystyle \frac{1}{x^{2n}+b^{2n}} = \sum_{k=0}^{2n-1} \frac{A_{k}}{x-be^{i \pi (2k+1)/(2n)}}$

where $\displaystyle A_{k} = \text{Res} \left[ \frac{1}{z^{2n}+b^{2n}}, be^{i \pi (2k+1) /(2n)} \right] = - \frac{1}{2nb^{2n-1}} e^{i \pi (2k+1)/(2n)}$

I'm going to denote $\displaystyle \frac{2k+1}{2n} \pi$ by $\displaystyle \lambda_{k}$

Then

$\displaystyle \frac{1}{x^{2n}+b^{2n}} = - \frac{1}{2nb^{2n-1}} \sum_{k=0}^{2n-1} \frac{e^{i \lambda_{k}}}{x-be^{i \lambda_{k}}}= - \frac{1}{2nb^{2n-1}} \sum_{k=0}^{n-1} \Big( \frac{e^{i \lambda_{k}}}{x-be^{i \lambda_{k}}} + \frac{e^{-i \lambda_{k}}}{x-be^{-i \lambda_{k}}} \Big)$

$\displaystyle = \frac{1}{nb^{2n-1}} \sum_{k=0}^{n-1} \frac{b- x \cos (\lambda_{k})}{x^{2}-2bx \cos( \lambda_{k})+b^{2}}$

$\displaystyle = \frac{1}{nb^{2n-1}} \sum_{k=0}^{n-1} \frac{b - x \cos (\lambda_{k})}{(x-b \cos \lambda_{k})^{2} + b^{2} \sin^{2} (\lambda_{k})}$

We also need the following two results which are fairly easy to derive using other known results.

$\displaystyle \int_{-\infty}^{\infty} \frac{\cos rx}{(x-b)^{2}+a^{2}} \ dx = \frac{\pi}{a} e^{-ar} \cos(br)$

$\displaystyle \int_{-\infty}^{\infty} \frac{x \cos (rx)}{(x-b)+a^{2}} \ dx = \frac{\pi}{a} e^{-ar} \left( b \cos br - a \sin br \right)$

Then

$\displaystyle \int_{-\infty}^{\infty} \frac{\cos (ax)}{x^{2n}+b^{2n}} \ dx = \frac{1}{nb^{2n-1}} \sum_{k=0}^{n-1} \left[ b \int_{-\infty}^{\infty} \frac{\cos (ax)}{(x-b \cos \lambda_{k})^{2}+ b^{2} \sin^{2} (\lambda_{k})} - \cos \lambda_{k} \int_{-\infty}^{\infty} \frac{ x \cos(ax)}{(x-b \cos \lambda_{k})^{2} + b^{2} \sin^{2} \lambda_{k}} \ dx \right]$

$\displaystyle = \frac{1}{nb^{2n-1}} \sum_{k=0}^{n-1} \Big[ \frac{\pi}{\sin \lambda_{k}} e^{-ab \sin \lambda_{k}} \cos(ab \cos \lambda_{k}) - \frac{\pi \cos \lambda_{k} }{b \sin \lambda_{k}} e^{-ab \sin \lambda_{k}} \Big( b \cos(\lambda_{k}) \cos( ab \cos \lambda_{k}) - b \sin (\lambda_{k}) \sin (ab \cos \lambda_{k} \Big) \Big]$

$\displaystyle = \frac{\pi}{nb^{2n-1}} \sum_{k=0}^{n-1} \frac{1}{\sin \lambda_{k}} e^{-ab \sin \lambda_{k}} \Big( \cos(ab \cos \lambda_{k}) - \cos(\lambda_{k}) \cos(ab \cos \lambda_{k}+\lambda_{k}) \Big)$

$\displaystyle = \frac{\pi}{nb^{2n-1}} \sum_{k=0}^{n-1} \frac{1}{\sin \lambda_{k}} e^{-ab \sin \lambda_{k}} \Big( \cos( ab \cos \lambda_{k}) - \cos^{2}(\lambda_{k}) \cos(ab \cos \lambda_{k}) + \cos (\lambda_{k}) \sin(\lambda_{k}) \sin(ab \cos \lambda_{k}) \Big)$

$\displaystyle = \frac{\pi}{nb^{2n-1}} \sum_{k=0}^{n-1} e^{-ab \sin \lambda_{k}} \Big( \sin(\lambda_{k}) \cos(ab \cos \lambda_{k}) + \cos(\lambda_{k}) \sin(ab \cos \lambda_{k}) \Big)$

$\displaystyle = \frac{\pi}{nb^{2n-1}} \sum_{k=0}^{n-1} e^{-ab \sin (\lambda_{k})} \sin \left(ab \cos \lambda_{k} + \lambda_{k} \right)$

Re: Integration Contest - Season 4

Sat Mar 22, 2014 4:42 pm
galactus
Global Moderator

Posts: 902
Very nice RV meister

I tried that one pretty much the same way. But man it was tedious

I will post my workings. It's quite a lot of typing though. Looks like you done a lot of typing too. Nice job.

Can I ask you something?. Have you ever done this one using contours?.

$\displaystyle \int_{0}^{\infty}\frac{\cos(x^{2})-\cos(x)}{x}dx$

Re: Integration Contest - Season 4

Sat Mar 22, 2014 9:00 pm
Random Variable Integration Guru

Posts: 381
Cody

I've never even considered trying to evaluate that integral using contour integration.

Can you get Maple or Mathematica to evaluate $\displaystyle \int_{0}^{\infty} \frac{\cos x}{1+x^{2n}}$ for a value of $n$ greater than or equal to $4$ and return something that doesn't look totally ridiculous?

Re: Integration Contest - Season 4

Sun Mar 23, 2014 11:18 am
galactus
Global Moderator

Posts: 902
Hey RV:

I done as you asked and ran that through Maple and Wolfram (I do not have regular Mathematica on my computer).

Wolfram returned nothing. (But, then again, it does not have the power as that of the regular computer-installed Mathematica).

Maple returned, as you said, gobbledy-gook in terms of the sum of the roots of some 2n degree polynomial and the Si and Ci functions. In other words, a mess. It may as well return nothing as far as a closed form goes. It returns numerical answers, but that's not what we want.

i.e. for n = 4, Maple gave $\displaystyle -\sum_{n=Root \;\ of \;\ z^{8}+1}\left[\frac{\pi \sin(n)}{16n^{7}}+\frac{Si(n)\sin(n)+Ci(-n)\cos(n)}{n^{7}}\right]$.

The larger n becomes, the ickier it gets.

This is what you worked out with a = 1 and b = 1 in the problem you just posted the solution to. See, it's cool when this happens. You solve a problem the fancy math engines can not.

I will try to take the time to post my workings for that one. It is lengthy. I tried contours as well.

Re: Integration Contest - Season 4

Sun Mar 23, 2014 12:37 pm
galactus
Global Moderator

Posts: 902
Well, RV, I will go ahead and post my workings to that one. All of the little intricacies were cumbersome. For me anyway. But, I got it assembled.

What I have here is for a = 1.

Write $\displaystyle \int_{C}\frac{e^{iz}}{b^{2n}+z^{2n}}dz$, where C is a typical semicircle in the UHP.

The poles are at $\displaystyle z=be^{\frac{(2k+1)\pi i}{2n}}$ and are equally spaced at $\displaystyle \frac{\pi}{n}$ intervals.

The poles that contribute are in the UHP. So, I am not going to bother writing out the others because they don't contribute anything anyway.

$\displaystyle \sum_{k=0}^{n-1}\frac{1}{2n\left(be^{\frac{(2k+1)}{2n}\pi i}\right)^{2n-1}\left(z-be^{\frac{(2k+1)\pi i}{2n}}\right)}$

The poles in the contour then contribute, due to the sum of the residues, in total (I hope this is correct. It appears to be. One itty-bitty negative or whatever and it's blown out of the water ):

$\displaystyle -\sum_{k=0}^{n-1}\frac{\pi i}{nb^{2n-1}}e^{-b\sin(\frac{2k+1}{2n}\pi}\left[\cos\left(\frac{2k+1}{2n}\pi +b\cos(\frac{2k+1}{2n}\pi)\right)+i\cdot \sin\left(\frac{2k+1}{2n}\pi +b\cos(\frac{2k+1}{2n}\pi)\right)\right]$.

The outside arc goes to 0 as R goes to infinity because of the $\displaystyle e^{-R\sin(t)}$ thing. So, I will not write it up.

Assembling:

$\displaystyle \int_{-\infty}^{0}\frac{e^{ix}}{b^{2n}+x^{2n}}dx+\int_{0}^{\infty}\frac{e^{ix}}{b^{2n}+x^{2n}}dx+\int_{0}^{\pi}\frac{e^{iRe^{it}}iRe^{it}}{b^{2n}+R^{2n}e^{2nit}}dt=$
$\displaystyle \sum_{k=0}^{n-1}\frac{\pi i}{nb^{2n-1}}e^{-b\sin(\frac{2k+1}{2n}\pi}\left[\cos\left(\frac{2k+1}{2n}\pi +b\cos(\frac{2k+1}{2n}\pi)\right)+i\cdot \sin\left(\frac{2k+1}{2n}\pi +b\cos(\frac{2k+1}{2n}\pi)\right)\right]$

Now, on the left, let $\displaystyle x\to -x$ in the first integral and add to the second one. Thus, giving:

$\displaystyle \int_{0}^{\infty}\frac{\cos(x)}{b^{2n}+x^{2n}}dx=\sum_{k=0}^{n-1}\frac{\pi}{2nb^{2n-1}}e^{-b\sin(\frac{2k+1}{2n}\pi)}\sin\left(\frac{2k+1}{2n}\pi +b\cos(\frac{2k+1}{2n}\pi)\right)$

There certainly could be a typo in all of that anyway.

For example, if n = 4, then the result is:

$\displaystyle \frac{\pi}{4}e^{-\sqrt{\sqrt{2}+2}\sin(\frac{\sqrt{2-\sqrt{2}}}{2}+\frac{3\pi}{8})/2}+\frac{\pi}{4}e^{-\sqrt{2-\sqrt{2}}\sin(\frac{\sqrt{2+\sqrt{2}}}{2}+\frac{\pi}{8})/2}$

which, I am sure, can be written in various forms, and appears to check numerically.

BTW, here is another general one to try with or without contours. Whatever you wish.

$\displaystyle \int_{0}^{\infty}\frac{\cos(x)}{(b^{2}+x^{2})^{n+1}}dx$

Maybe a little tougher?. Anyway, the next time an integral like this pops up on SE or wherever, one can just show the general case and maybe show off a little bit :):).

Previous


text:
## Algebra and Trigonometry 10th Edition

$x \gt 0\\y \gt x \\ \pi y^2 \gt \pi x^2+10$.
The radius of the smaller circle is positive; this means that $x \gt 0$ The radius $y$ is for the larger circle; this means that $y \gt x$ The area of the larger circle is at least $10 \ sq \ unit$ greater; this means that $\pi y^2 \gt \pi x^2+10$. Thus, the system is: $x \gt 0\\y \gt x \\ \pi y^2 \gt \pi x^2+10$.


text:
# Commutators in 2 dimensions

1. Nov 29, 2007

### Sunshine

1. The problem statement, all variables and given/known data

Simplify the commutator [A,B] and give the expectation value of [A,B] in the ground state for an isotropic harmonic oscillator (mass m) that has the energy $$\hbar \omega /2$$ when
$$A = xp_x$$
$$B = y$$

2. Relevant equations
$$[AB,C] = A[B,C] + [A,C]B$$
$$[p_i,x_j] = i\hbar\delta_{ij}$$

3. The attempt at a solution
$$[xp_x,y] = x[p_x,y] + [x,y]p_x$$ (first relation)
$$x[p_x,y] = 0$$ (second relation)
Last term with test function f(x)
$$[x,y]p_x f(x) = xy\frac\hbar i \dfrac{\partial}{\partial x}f(x) - yx \frac\hbar i \dfrac{\partial}{\partial x}f(x) = 0 ?$$

I have a feeling that 0 isn't the answer, since I have to find the expectation value as well. If the last equation doesn't become 0 but the middle equation is the most simplified answer, I don't know how to find an expectation value that isn't equal to 0 (because I get that when I put it into the usual integral for expectation value)

2. Nov 29, 2007

### nrqed

the commutator of the A and B you gave above is indeed zero which is obvious because y commutes with p_x and with y as well.

Are you sure that A is not x p_y instead of x p_x?? Because then the calculation would be more interesting.

3. Nov 29, 2007

### Sunshine

yes, it is $$xp_x$$.

So the expectation value is 0 as well?

4. Nov 29, 2007

### nrqed

Yes, it is zero.

Strange question :uhh:

5. Nov 29, 2007

### Sunshine

I'm confused because it seems too easy to count expectation value for 0. Especially when question 2 (out of 3 in the assignment) has a trivial expectation value too. Doing it for 3rd time just to make sure I'm not missing something. If it becomes zero again, I'll be almost certain that something is wrong.

6. Nov 29, 2007

### nrqed

is that a question from a book or was it typed up by your professor? I would think that he/she may have made a mistake or that there is a typo in the book.
(or do they want the expectation value of AB instead of [A,B]?)

7. Nov 29, 2007

### Sunshine

The professor took the question from an unknown book, typed it into word and emailed it. I don't think it's a typo, and I'd never dare to ask if it was. He wouldn't admit his mistake anyway ;)

ok, second question is like the first one with a few differences:

Now I have
$$A = L_x$$
$$B = z_x^2$$

And it's not a harmonic oscillator, but a hydrogen atom in ground state. My solution:

$$[L_x, zx^2] = [yp_z - zp_y, zx^2]$$
$$=[yp_z,zx^2]-[zp_y,zx^2]$$

First term: $$y(\mathbf{[p_z,z]x^2} + z[p_z,x^2])+([y,z]x^2 + z[y,x^2])p_z$$
Second term:$$z([p_y,z]x^2 + z[p_y,x^2])+([z,z]x^2 + z[z,x^2])p_y$$

Only the term in bold is non-zero, which gives:

$$y[p_z,z]x^2 = -i\hbar y x^2$$

The squared absolute value of wavefunction for hydrogen atom in ground state:
$$|\Psi |^2 = \dfrac{1}{\pi a^3} e^{-2r/a}$$

Hence, the expectation value is (changing to spherical coordinates):
$$-ih<yx^2> = ih<r^2 \sin^3\theta \cos^2\phi \sin\phi > = C \int_0^{2\pi} \int_0^\pi \int_0^a (r^2 \sin^3\theta \cos^2\phi \sin\phi) e^{-2r/a} r^2 \sin \theta dr d\theta d\phi$$

This is 0! (Tried with Mathematica) I'm not sure about the interval for r; if it's from 0 to a or from 0 to infinity. In the last case the integral doesn't converge.

Can you see where I'm going wrong?

Last edited: Nov 29, 2007
8. Nov 30, 2007

### nrqed

It (almost) all looks right.

Within doing any integral it is clear that the integration will give zero. The ground state wavefunction is spherically symmetric so it's even under $$x,y,z \rightarrow -x,-y,-z$$ but the integrand xy^2 is obviously odd, so the integral must vanish.

The only wrong thing I can see is that you should be able to integrate r from 0 to infinity (the exponential will make any integral with powers of r converge). But that does not change the final conclusion, the integral gives zero.

9. Nov 30, 2007

### Sunshine

You are right, the last integral becomes 0 as well -I accidentally made the assumption that a<0 istead of a>0 in Mathematica last time.

So the conclusion is that the expectation value of both commutators is 0?

10. Nov 30, 2007

### nrqed

I am not sure what you mean by last integral but if you mean the integral over r, then no, it does not give zero by itself. And I am not sure what "a" is. You mean the Bohr radius? No, the radial integral should be from 0 to infinity, period. But the angular integral will give zero so that the final result is zero.

Looks like it to me.


text:
anonymous 5 years ago solve (two answers) (x-3)^2/3=9

$(x-3)=9^{\frac{3}{2}}=27$ $x=30$ or $x-3=-27$ $x=-24$


text:
Fundamentals of Physics Extended (10th Edition)

$F = 6. 15$ $N$
The force exerted on the top link can be calculated by adding the weights of all links to the force required to accelerate all the links at $2.50m/s^{2}$, $F=(m_{1}+m_{2}+m_{3}+m_{4}+m_{5})g + (m_{1}+m_{2}+m_{3}+m_{4}+m_{5})a$ $F=(m_{1}+m_{2}+m_{3}+m_{4}+m_{5} ) (g +a)$ $F =[(0.500kg) (9.8m/s^{2}+2.5m/s^{2})]=6.15$ $N$


text:
Deere & Co., the Moline, IL, agricultural equipment giant, has been using rapid prototyping (RP) for about a decade, and is now finishing a wider examination of its entire product development process. Deere is seeking to cut development time by as much as 50 percent, and is looking at how the different prototyping techniques can help that along. Using solid modeling software, as Deere did, is critical for engineering teams wishing to use RP technology. Deere has expended considerable effort recently developing a new product-delivery process (PDP). Such processes can be difficult to hammer out and then follow faithfully, especially in today’s work environment of lean staffs and greater profit pressure. Although RP materials have improved significantly since the first stereolithography machine came into use a little over 10 years ago, none of the current materials to date have the durability required for use on John Deere products sold to customers.

## Article

Throughout the 1990s, engineering departments and divisions worldwide have been acclimating and refining the processes of rapid prototyping. New methods have come to the fore, and the cost of the original methods has come down dramatically. Systems have even made the proverbial march to the desktop.

Deere & Co., the Moline, Ill., agricultural equipment giant, has been using RP for about a decade, and is now finishing a wider examination of its entire product development process. Deere is seeking to cut development time by as much as 50 percent, and is looking at how the different prototyping techniques can help that along.

We talked with John K. Lawson, Deere's senior VP for engineering, information, and technology, about RP and the engineering process.

Lawson reports that rapid prototyping of physical models has been important in Deere's product design cycle, but today the role of virtual prototyping—that is, computer modeling-is growing faster, and the company is finding advantages in working with computer models at ever-earlier stages of design.

Senior VP John K. lawson of Deere: "We move from solving problems to preventing them."

## Proving Rapid Prototyping's Worth

John Deere has been using what is commonly referred to as rapid prototyping to produce physical models since 1990. At that time, the new technology was seen as another potentially valuable tool to improve the development of John Deere products, Lawson said, but its worth had to be proved. Deere executives recognized that RP, like many emerging technologies, would evolve and develop over time, that some of the processes would mature into a lasting, robust technology, and others would disappear. With this in mind, Deere decided not to in- vest in the ownership of RP equipment, but opted instead to purchase RP models and replicas from outside service bureaus.

Deere appointed an RP facilitator to introduce the new technology at various product development sites, to assist specific design teams in rapid prototyping, and to conduct in-house seminars. The approach enabled Deere to involve about 200 designers, engineers, and other design-team members in the use of RP during the early 1990s.

Since Deere had decided not to purchase RP equipment, a sm60th transition from in-house engineering to the service bureaus was important. Integrating a computer modeling program-in this case, Parametric Technology Corp.'s Pro/Engineer-into the design process made it easier for designers to order RP models from the company's outside supplier.

Using solid modeling software, as Deere did, is critical for engineering teams wishing to use RP technology, according to Terry Wohlers, an industry consultant and analyst based in Fort Collins, Colo. Nearly all RP systems start from a stereolithography (STL) file, which is the de facto standard format for RP systems and is supported by almost all CAD programs. These STL files must be what Wohlers terms "watertight," so the RP system does not choke on holes and gaps in the data.

Lawson said that only in the last year or so has John Deere purchased any equipment for in-house building of physical RP models. Now, one Genisys desktop machine from Stratasys and one Actua 2100 desktop modeler from 3D Systems are in use at Deere's Product Engineering Center in Waterloo, Iowa, where tractor and engine design teams are making concept parts. Another Genisys machine, at the Deere Technical Center in Moline, has been used as a more centralized source by different product-design centers.

Deere engineers exploited their CAE software, including PTC's Pro/Engineer, to virtually prototype the components of this tractor control panel.

## Establshing Engineering Procedure

Deere has expended considerable effort recently developing a new product-delivery process (PDP), Lawson said. Such processes can be difficult to hammer out and then follow faithfully, especially in today's work environment of lean staffs and greater profit pressure. But for the same reasons, a clear process that everyone follows is vital.

The Deere process is "a structure that defines the steps and approvals that are necessary to develop a new product," Lawson said. But the process does not enumerate how many engineers or how much time each step can take. Products provided by the four Deere equipment divisions vary considerably in their complexity and size, he said, so for this and other reasons there can't be specific time frames or manpower limits for product development.

"In a very large complex product PDP, we may have as many as 75 engineers over a period of three to four years," he said. "A smaller program may have 10. engineers and be complete within nine months." In the end, Deere wants to examine the product development cycle and then shorten it by 30. to 50. percent. Of course, Lawson said, the reduction in time has to occur along with improvements in product and process quality.

The first steps of the product delivery process break out in this order: strategic business planning, in which the customer and market are the first consideration (an engineering maxim); program definition, where potential design and manufacturing concepts are laid out; concept evaluation and program specification; and program development, when detailed design and development are finished.

In the fifth stage, confirmation and implementation, the manufacturing prototypes are made on the assembly line. The product delivery process is in a continuous sixth-and final-stage during manufacturing and delivery, when feedback is gathered for continuous improvement.

It is the fourth stage, program development, where rapid prototyping often falls, though Wohlers, the industry consultant, argued that engineering time is reduced an order of magnitude for each stage earlier in the game that RP is brought in.

Deere has found that earlier product design pays dividends, but not with rapid prototyping as much as with virtual prototyping, which is made possible by the latest CAE software.

According to Lawson, the development work with "hard prototypes" is typically done in the fifth phase, confirmation and implementation. "If we develop complete virtual prototypes in phases two and three," he said, "we move from solving problems to preventing problems."

However, brand-new technology requires a second protocol. "We are in the process of developing an enterprise tec hnology delivery process," Lawson said. "We have found that developing new technology must be done outside of a product PDP. And that TDP must be complete before inclusion in a PDP."

The technology delivery planat Deere includes technology identification, research, assessment, and development. The concept is to use only "proven" technology in product programs.

Components of this John Deere 8000 series tractor panel were rapid prototyped during the design phase.

## Building Fewer Prototypes

In some instances, advances in software appear to be presenting a serious challenge to the rationale for rapid prototyping. "At Deere, products are often quite large and involve various metal castings, sheet metal welded components, and plastic parts," Lawson said, "so what is commo nly referred to as RP-stereolithography, selective laser sintering, laminated object manufacturing, and soon-can only be applied to a portion of our designs."

As the product delivery protocol indicates, "We are moving toward electronic-design iterations, using product visualization software and virtual reality," he said, "and we are building fewer physical prototypes."

Whether a particular part is suitable for RP, or needs to be prototyped at all, is a decision left to members of the specific design team. "With Pro/Engineer as the foundation for our designs, it really makes no difference to the designer," Lawson said. Again, that 's beca use the Prol E files are easily transferred into RP-suitable STL files. The cost model is encouraging the drive away from physical prototypes and toward virtual prototyping, also. Most common full-size, robust RP machines, fully installed, cost $300,000 to$500,000 each, Lawson said.

"And in order to match specific part needs with the most suited technology, you might have to spend over a million dollars, in addition to the cost of a climate-controlled room, auxiliary equipment, and so on," he said. The number of physical models being built each year, even at a large company like Deere, "would not justify such expenditure."

## 'Desktop-Like' Systems

So, Deere employs a few of the new" desktop-like" systems and uses the mature process in place for sending CAD files out to an RP service bureau when that is deemed necessary. But th e emphasis is on virtual prototyping.

Lawson brought to light one reason that rapid prototyping has reached this more mature stage. RP still cannot be small-lot computer-aided manufac turing, because "the RP equipment does not yet produce a part with true engineering material properties and with the speed and cost near that of CNC machining, injection molding, or other manufacturing processes," Lawson said.

Although RP materials have improved significantly since the first stereolithography machine came into use a little over 10 years ago, none of the current materials to date have the durability required for use on John Deere products sold to customers," he added. In several isolated cases, RP components have been used for short-term field testing.

There are a number of ways that the field of rapid prototyping can broaden. Expect a greater number of options th anever, especially in software. And, as Wohlers has reported, expect to hear that engineers from service bureaus are expanding their capabilities to be sure they can continue to service the greater engineering community.


text:
Indefinite Integrals

# Indefinite Integrals

We've look at The Fundamental Theorem of Calculus Part 1 and The Fundamental Theorem of Calculus Part 2 to calculate definite integrals. There are some cases in which we don't want to calculate an integral on an interval $[a, b]$, but instead, we may be more interested with the general integral itself.

 Definition: If $f$ is a function, then the Indefinite Integral of $f$ denoted $\int f(x) \: dx = F(x) + C$ where $F(x) + C$ is any antiderivative of $f$.

We should note that indefinite integration is analogous to that of anti-differentiation, and as a result, we obtain the following list of indefinite integrals:

• If $k$ is a real number, then $\int k f(x) \: dx = k \int f(x) \: dx$.
• If $k$ is a real number, then $\int k \: dx = kx + C$.
• If $n ≠ 1$, then $\int x^n \: dx = \frac{x^{n+1}}{n + 1} + C$.
• $\int \sin x \: dx = -\cos x + C$.
• $\int \cos x \: dx = \sin x + C$.
• $\int \sec ^2 x \: dx = \tan x + C$.
• $\int \sec x \tan x \: dx = \sec x + C$.
• $\int \csc x \cot x \: dx = -\csc x + C$.
• $\int a^x \: dx = \frac{a^x}{\ln a} + C$.
• $\int e^x \: dx = e^x + C$.
• $\int \frac{1}{x} \: dx = \ln \mid x \mid + C$.
• $\int \frac{1}{\sqrt{1 - x^2}} \: dx = \sin ^{-1} x + C$.
• $\int \frac{1}{x^2 + 1} \: dx = \tan ^{-1} x + C$.
 Remark: Like when dealing with antiderivatives, it is important to always add a $+C$ for indefinite integrals! Many courses will deduct marks if this detail is forgotten.