Hey guys,
I am doing a program design of a porous medium magnetic fluid model. Now there is a problem that needs help. The following is my program code。
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
#setting parameters
M = 1.0
Lambda = 1.0
K = 1.0 # K^{-1} in paper
niu_w = 1.0
niu_f = 2.0
miu = 1.0
tau = 1.0
beta = 1.0
B = 1.0 # SAV constant for ensuring sqrt{} meaningful
k_w = 0
k_f = 1.0
EPS = 0.5
# Artifical force
f_phi = Expression('cos(pi*x[0])*cos(pi*x[1])*cos(t) + 2*M*pow(pi, 2)*cos(pi*x[0])*cos(pi*x[1])*sin(t) - pi*pow(cos(pi*x[0]), 2)*sin(pi*x[1])*pow(sin(t), 2)*cos(pi*(x[1] + 0.5)) - pi*cos(pi*x[1])*pow(sin(pi*x[0]), 2)*pow(sin(t), 2)*sin(pi*(x[1] + 0.5))', degree=2, M=M, t=0)
f_w = Expression('cos(pi*x[0])*cos(pi*x[1])*sin(t) - 2*Lambda*pow(pi, 2)*cos(pi*x[0])*cos(pi*x[1])*sin(t) - (Lambda*cos(pi*x[0])*cos(pi*x[1])*sin(t)*(pow(cos(pi*x[0]),2)*pow(cos(pi*x[1]),2)*pow(sin(t),2)- 1))/pow(EPS,2)', degree=2, Lambda=Lambda, EPS=EPS, t=0)
f_u = Expression(('2*sin(t)*(2*x[1]-1)-miu*sin(t+x[0])*sin(t)+sin(pi*x[0])*cos(t)*sin(pi*(x[1]+0.5))-pi*cos(pi*x[0])*pow(cos(pi*x[1]),2)*sin(pi*x[0])*pow(sin(t),2)+K*sin(pi*x[0])*sin(t)*sin(pi*(x[1]+0.5))*(1/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS)+1)+1)',
'2*sin(t)*(2*x[0]-1)-miu*sin(t+x[1])*sin(t)+cos(pi*x[0])*cos(t)*cos(pi*(x[1]+0.5))-pi*pow(cos(pi*x[0]),2)*cos(pi*x[1])*sin(pi*x[1])*pow(sin(t),2)+K*cos(pi*x[0])*sin(t)*cos(pi*(x[1]+0.5))*(1/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS)+1)+1)'), degree=2, K=K, EPS=EPS, miu=miu, t=0)
f_m = Expression(('cos(t + x[1]) + (sin(t + x[1]) - (x[1]*sin(t))/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS) + 1))/tau - beta*sin(t + x[0])*(x[1]*sin(t + x[0])*sin(t) - sin(t + x[1])*sin(t)*(x[0] - 0.5)) + cos(pi*x[0])*cos(t + x[1])*sin(t)*cos(pi*(x[1] + 0.5))',
'cos(t + x[0]) + (sin(t + x[0]) - (sin(t)*(x[0] - 0.5))/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS) + 1))/tau + beta*sin(t + x[1])*(x[1]*sin(t + x[0])*sin(t) - sin(t + x[1])*sin(t)*(x[0] - 0.5)) + sin(pi*x[0])*cos(t + x[0])*sin(t)*sin(pi*(x[1] + 0.5))'), degree=2, EPS=EPS, beta=beta, tau=tau, t=0)
h_a = Expression(('sin(t+x[1])+x[1]*sin(t)', 'sin(t+x[0])+(x[0]-0.5)*sin(t)'), degree=2, t=0)
h_b = Expression(('cos(t+x[1])+x[1]*cos(t)', 'cos(t+x[0])+(x[0]-0.5)*cos(t)'), degree=2, t=0)
# Exact solution
phi_e = Expression('sin(t)*cos(pi*x[0])*cos(pi*x[1])', degree=2, t=0) #exact solution of theta
w_e = Expression('sin(t)*cos(pi*x[0])*cos(pi*x[1])', degree=2, t=0) #exact solution of omega
u_e = Expression(('sin(t)*sin(pi*x[0])*sin(pi*(x[1]+0.5))', 'sin(t)*cos(pi*x[0])*cos(pi*(x[1]+0.5))'), degree=2, t=0) # exact solution of u
p_e = Expression('sin(t)*(2*x[0]-1)*(2*x[1]-1)', degree=2, t=0) #exact solution of p
m_e = Expression(('sin(t+x[1])', 'sin(t+x[0])'), degree=2, t=0) #exact solution of m
varphi_e = Expression('(x[0]-0.5)*x[1]*sin(t)', degree=2, t=0) #exact solution of varphi
h_e = Expression(('x[1]*sin(t)', '(x[0]-0.5)*sin(t)'), degree=2, t=0) #exact solution of h
# Denfine dt and mesh
T = 0.5
m = 32 #time step
dt = 1.0/m
k = Constant(dt) # \delta t
n = 32 # space mesh
mesh = UnitSquareMesh(n, n)
domain = 1.0 * 1.0 #domain
# Denfine the finite element and the mixed functionspace
P1 = FiniteElement("CG", mesh.ufl_cell(), 1) # P1 element
TH1 = MixedElement([P1, P1]) # Mixed element P1*P1
W1 = FunctionSpace(mesh, TH1) # Mixed functionspace P1*P1
# Denfin the functionspace
P1Space = FunctionSpace(mesh, "CG", 1) #P1 Space
P1_dgSpace = VectorFunctionSpace(mesh, "DG", 1) #P1-dg Space
P1VectorSpace = VectorFunctionSpace(mesh, "CG", 1) #P1-vector Space
P2Space = FunctionSpace(mesh, "CG", 2) #P2 Space
# solve the average
int = Function(P2Space)
Unit = Constant(1.0)
average = (1.0/domain)*dot(int, Unit)*dx #solve (1/|Omega|)*int_{Omega} ... dx
#Define function
def F(a):
return (1.0/4)*(pow(a, 2)-1)*(pow(a, 2)-1) # define the F(phi)
def f(a):
return a*(pow(a, 2)-1) #define the f(fai)
def vcross(x, y):
return x[0] * y[1] - x[1] * y[0] # 向量叉乘
# def h(a):
# return grad(a) #define the h=grad(fai)
def niu(a):
return niu_w + (niu_f-niu_w)*(1.0/(1.0+exp(-(a/EPS)))) #define the \nu(phi)
def kappa(a):
return k_w+(k_f-k_w)*(1.0/(1.0+exp(-(a/EPS)))) #define the \kappa(phi)
def intergal(a):
return dot(F(a), 1)*dx #define int_{Omega} F(phi)*dx
# Define the initial value
phi0 = interpolate(phi_e, P1Space) # phi0
w0 = interpolate(w_e, P1Space)
# w0_trial = TrialFunction(P1Space)
# y0_test = TestFunction(P1Space)
# w0 = Function(P1Space) # w0
# aw0 = dot(w0_trial, y0_test)*dx
# Lw0 = Lambda*dot(grad(phi0), grad(y0_test))*dx + (Lambda/(EPS*EPS))*dot(f(phi0), y0_test)*dx
# solve(aw0 == Lw0, w0) # solve new w0
u0 = interpolate(u_e, P1_dgSpace) # u0
p0 = Constant(0) # p0
m0 = interpolate(m_e, P1VectorSpace) # m0
ha0 = interpolate(h_a, P1VectorSpace)
varphi0_trial = TrialFunction(P2Space)
psi0_test = TestFunction(P2Space)
varphi0 = Function(P2Space)
avarphi0 = dot(grad(varphi0_trial), grad(psi0_test))*dx
Lvarphi0 = -dot(m0, grad(psi0_test))*dx + dot(ha0, grad(psi0_test))*dx
solve(avarphi0 == Lvarphi0, varphi0)
int.assign(varphi0)
value0 = assemble(average)
varphi0.vector()[:] = varphi0.vector() - value0
r0 = sqrt( (1.0/(EPS*EPS))*assemble(intergal(phi0)) + B ) # r0
print('r0= ', r0)
# Define the old value
phi_old = Function(P1Space) # phi^{n}
w_old = Function(P1Space) # w^{n}
u_old = Function(P1_dgSpace) # u^{n}
p_old = Function(P2Space) # p^{n}
m_old = Function(P1VectorSpace) # m^{n}
varphi_old = Function(P2Space) # varphi^{n}
r_old = 1.0
phi_a = Function(P1Space) # phi_a^{n+1}
phi_b = Function(P1Space) # phi_B^{n+1}
w_a = Function(P1Space) # w_a^{n+1}
w_b = Function(P1Space) # w_b^{n+1}
utilde_a = Function(P1_dgSpace) # save utilde_a
utilde_b = Function(P1_dgSpace) # save utilde_b
varphi_a = Function(P2Space) #save varphi_a
varphi_b = Function(P2Space) #save varphi_b
m_a = Function(P1VectorSpace) #save m_a
m_b = Function(P1VectorSpace) #save m_b
phi = Function(P1Space) # save phi^{n+1}
w = Function(P1Space) # save w^{n+1}
utilde = Function(P1_dgSpace) #save utilde^{n+1}
varphi = Function(P2Space) #save varphi^{n+1}
m = Function(P1VectorSpace) #save m^{n+1}
p = Function(P2Space) #save p^{n+1}
u = Function(P1_dgSpace) #save u^{n+1}
S_n = 1.0
r = 1.0
# Define trial and test function
(phia_trial, wa_trial) = TrialFunctions(W1)
(X_test, Y_test) = TestFunctions(W1)
G1 = Function(W1) # save phi_a^{n+1} and w_a^{n+1}
aphi_a = (1.0/k)*dot(phia_trial, X_test)*dx + M*dot(grad(wa_trial), grad(X_test))*dx + dot(wa_trial, Y_test)*dx - Lambda*dot(grad(phia_trial), grad(Y_test))*dx
Lphi_a = (1.0/k)*dot(phi_old, X_test)*dx + dot(f_phi, X_test)*dx + dot(f_w, Y_test)*dx
(phib_trial, wb_trial) = TrialFunctions(W1)
G2 = Function(W1) # save phi_b^{n+1} and w_b^{n+1}
aphi_b = (1.0/k)*dot(phib_trial, X_test)*dx + M*dot(grad(wb_trial), grad(X_test))*dx + dot(wb_trial, Y_test)*dx - Lambda*dot(grad(phib_trial), grad(Y_test))*dx
Lphi_b = (1.0/S_n)*dot(u_old*phi_old, grad(X_test))*dx + (1.0/S_n)*(Lambda/(EPS*EPS))*dot(f(phi_old), Y_test)*dx
Phn_trial = TrialFunction(P1_dgSpace)
n_test = TestFunction(P1_dgSpace)
Ph_n = Function(P1_dgSpace)
aPh = dot(Phn_trial, n_test)*dx
LPh = dot(grad(varphi_old), n_test)*dx
utildea_trial = TrialFunction(P1_dgSpace)
v_test = TestFunction(P1_dgSpace)
autilde_a = (1.0/k)*dot(utildea_trial, v_test)*dx + K*dot(niu(phi_old)*utildea_trial, v_test)*dx
Lutilde_a = (1.0/k)*dot(u_old, v_test)*dx - dot(grad(p_old), v_test)*dx + dot(f_u, v_test)*dx
utildeb_trial = TrialFunction(P1_dgSpace)
autilde_b = (1.0/k)*dot(utildeb_trial, v_test)*dx + K*dot(niu(phi_old)*utildeb_trial, v_test)*dx
Lutilde_b = -(1.0/S_n)*dot(phi_old*grad(w_old), v_test)*dx + (1.0/S_n)*miu*dot(dot(m_old, nabla_grad(Ph_n)), v_test)*dx
varphia_trial = TrialFunction(P2Space)
psi_test = TestFunction(P2Space)
avarphi_a = (1.0/k)*dot(grad(varphia_trial), grad(psi_test))*dx + (1.0/tau)*dot(grad(varphia_trial), grad(psi_test))*dx + (1.0/tau)*dot(kappa(phi_old)*grad(varphia_trial), grad(psi_test))*dx + beta*dot(vcross(m_old, grad(varphia_trial)), vcross(m_old, grad(psi_test)))*dx
Lvarphi_a = (1.0/k)*dot(grad(varphi_old), grad(psi_test))*dx + (1.0/tau)*dot(h_a, grad(psi_test))*dx + dot(h_b, grad(psi_test))*dx - dot(f_m, grad(psi_test))*dx
varphib_trial = TrialFunction(P2Space)
avarphi_b = (1.0/k)*dot(grad(varphib_trial), grad(psi_test))*dx + (1.0/tau)*dot(grad(varphib_trial), grad(psi_test))*dx + (1.0/tau)*dot(kappa(phi_old)*grad(varphib_trial), grad(psi_test))*dx + beta*dot(vcross(m_old, grad(varphib_trial)), vcross(m_old, grad(psi_test)))*dx
Lvarphi_b = (1.0/S_n)*dot(dot(u_old, nabla_grad(m_old)), grad(psi_test))*dx
ma_trial = TrialFunction(P1VectorSpace)
z_test = TestFunction(P1VectorSpace)
am_a = (1.0/k)*dot(ma_trial, z_test)*dx + (1.0/tau)*dot(ma_trial, z_test)*dx
Lm_a = (1.0/k)*dot(m_old, z_test)*dx + (1.0/tau)*dot(kappa(phi_old)*grad(varphi_a), z_test)*dx + dot(f_m, z_test)*dx
mb_trial = TrialFunction(P1VectorSpace)
am_b = (1.0/k)*dot(mb_trial, z_test)*dx + (1.0/tau)*dot(mb_trial, z_test)*dx
Lm_b = (1.0/tau)*dot(kappa(phi_old)*grad(varphi_b), z_test)*dx - (1.0/S_n)*dot(dot(u_old, nabla_grad(m_old)), z_test)*dx + (1.0/S_n)*beta*dot(vcross(m_old, grad(varphi_old)), vcross(m_old, z_test))*dx
eta_a = (1.0/(2*S_n*pow(EPS, 2)))*dot(f(phi_old), phi_a)*dx - (k/(2*Lambda*S_n))*dot(u_old*phi_old, grad(w_a))*dx + (k/(2*Lambda*S_n))*dot(phi_old*grad(w_old), utilde_a)*dx \
- (k/(2*Lambda*S_n))*miu*dot(dot(m_old, nabla_grad(Ph_n)), utilde_a)*dx - (k/(2*Lambda*S_n))*miu*dot(dot(u_old, nabla_grad(m_old)), grad(varphi_a))*dx \
+ (k/(2*Lambda*S_n))*(miu/k_f)*dot(dot(u_old, nabla_grad(m_old)), m_a)*dx - (k/(2*Lambda*S_n))*(miu/k_f)*beta*dot(vcross(m_old, grad(varphi_old)), vcross(m_old, m_a))*dx
Eta_a = assemble(eta_a)
eta_b = (1.0/(2*S_n*pow(EPS, 2)))*dot(f(phi_old), phi_b)*dx - (k/(2*Lambda*S_n))*dot(u_old*phi_old, grad(w_b))*dx + (k/(2*Lambda*S_n))*dot(phi_old*grad(w_old), utilde_b)*dx \
- (k/(2*Lambda*S_n))*miu*dot(dot(m_old, nabla_grad(Ph_n)), utilde_b)*dx - (k/(2*Lambda*S_n))*miu*dot(dot(u_old, nabla_grad(m_old)), grad(varphi_b))*dx \
+ (k/(2*Lambda*S_n))*(miu/k_f)*dot(dot(u_old, nabla_grad(m_old)), m_b)*dx - (k/(2*Lambda*S_n))*(miu/k_f)*beta*dot(vcross(m_old, grad(varphi_old)), vcross(m_old, m_b))*dx
Eta_b = assemble(eta_b)
eta_c = (1.0/(2*S_n*pow(EPS, 2)))*dot(f(phi_old), phi_old)*dx
Eta_c = assemble(eta_c)
phi_trial = TrialFunction(P1Space)
x_test = TestFunction(P1Space)
aphi = dot(phi_trial, x_test)*dx
Lphi = dot(phi_a + r*phi_b, x_test)*dx
w_trial = TrialFunction(P1Space)
y_test = TestFunction(P1Space)
aw = dot(w_trial, y_test)*dx
Lw = dot(w_a + r*w_b, y_test)*dx
utilde_trial = TrialFunction(P1_dgSpace)
autilde = dot(utilde_trial, v_test)*dx
Lutilde = dot(utilde_a + r*utilde_b, v_test)*dx
varphi_trial = TrialFunction(P2Space)
avarphi = dot(varphi_trial, psi_test)*dx
Lvarphi = dot(varphi_a + r*varphi_b, psi_test)*dx
m_trial = TrialFunction(P1VectorSpace)
am = dot(m_trial, z_test)*dx
Lm = dot(m_a + r*m_b, z_test)*dx
p_trial = TrialFunction(P2Space)
q_test = TestFunction(P2Space)
ap = dot(grad(p_trial), grad(q_test))*dx
Lp = (1.0/k)*dot(utilde, grad(q_test))*dx + dot(grad(p_old), grad(q_test))*dx
u_trial = TrialFunction(P1_dgSpace)
au = dot(u_trial, v_test)*dx
Lu = dot(utilde, v_test)*dx - k*dot(grad(p), v_test)*dx + k*dot(grad(p_old), v_test)*dx
phi_old.assign(phi0)
w_old.assign(w0)
u_old.assign(u0)
p_old.assign(p0)
varphi_old.assign(varphi0)
m_old.assign(m0)
r_old = r0
# Work in cycle
t = dt
while t < T + 1.0e-8:
print('t= ', t)
f_phi.t = t
f_w.t = t
f_u.t = t
f_m.t = t
h_a.t = t
h_b.t = t
S_n = sqrt( (1.0/(EPS*EPS))*assemble(intergal(phi_old)) + B) # define S_n
solve(aphi_a == Lphi_a, G1) # solve phi_a^{n+1} and w_a^{n+1}
(phia_up, wa_up) = G1.split(True)
phi_a.assign(phia_up)
w_a.assign(wa_up)
solve(aphi_b == Lphi_b, G2) # solve phi_b^{n+1} and w_b^{n+1}
(phib_up, wb_up) = G2.split(True)
phi_b.assign(phib_up)
w_b.assign(wb_up)
solve(aPh == LPh, Ph_n) # solve the L^{2} projection P(h)
solve(autilde_a == Lutilde_a, utilde_a) # solve utilde_a^{n+1}
solve(autilde_b == Lutilde_b, utilde_b) # solve utilde_b^{n+1}
solve(avarphi_a == Lvarphi_a, varphi_a) # solve varphi_a^{n+1}
int.assign(varphi_a)
value1 = assemble(average)
varphi_a.vector()[:] = varphi_a.vector() - value1 # # varphi_a在[0,1]上的积分为0
solve(avarphi_b == Lvarphi_b, varphi_b) # solve varphi_b^{n+1}
int.assign(varphi_b)
value2 = assemble(average)
varphi_b.vector()[:] = varphi_b.vector() - value2 # varphi_b在[0,1]上的积分为0
solve(am_a == Lm_a, m_a) # solve m_a^{n+1}
solve(am_b == Lm_b, m_b) # solve m_b^{n+1}
Eta_a = assemble(eta_a)
Eta_b = assemble(eta_b)
Eta_c = assemble(eta_c)
r = (Eta_a - Eta_c + r_old) / (1.0 - Eta_b) #solve r^{n+1}
solve(aphi == Lphi, phi) # combine phi
solve(aw == Lw, w) # combine w
solve(autilde == Lutilde, utilde) # combine utilde
solve(avarphi == Lvarphi, varphi) # combine varphi
int.assign(varphi)
value3 = assemble(average)
varphi.vector()[:] = varphi.vector() - value3 # varphi^{n+1}在[0,1]上的积分为0
solve(am == Lm, m) # combine m
solve(ap == Lp, p) # solve p^{n+1}
int.assign(p)
value4 = assemble(average)
p.vector()[:] = p.vector() - value4 ## p^{n+1}在[0,1]上的积分为0
solve(au == Lu, u) # solve u^{n+1}
print('r=', r)
print('S_n=', S_n)
phi_old.assign(phi)
w_old.assign(w)
u_old.assign(u)
p_old.assign(p)
varphi_old.assign(varphi)
m_old.assign(m)
r_old = r
t = t + dt
# Compute the error
phi_e.t = T
u_e.t = T
p_e.t = T
m_e.t = T
varphi_e.t = T
h_e.t = T
error_phi_L2 = errornorm(phi_e, phi, 'L2')
error_phi_H1 = errornorm(phi_e, phi, 'H1')
error_uL2 = errornorm(u_e, u, 'L2')
error_pL2 = errornorm(p_e, p, 'L2')
error_pH1 = errornorm(p_e, p, 'H1')
error_mL2 = errornorm(m_e, m, 'L2')
h = grad(varphi)
error_hL2 = (h_e - h)**2*dx
E_h = sqrt(abs(assemble(error_hL2)))
phie = interpolate(phi_e, P1Space)
error_r = abs(r - sqrt((1.0/(EPS*EPS))*assemble(intergal(phie)) + B))
print('error_phi_L2 =', error_phi_L2)
print('error_phi_H1 =', error_phi_H1)
print('error_uL2 =', error_uL2)
print('error_pL2 =', error_pL2)
print('error_pH1 = ', error_pH1)
print('error_mL2 =', error_mL2)
print('error_hL2 = ', E_h)
print('error_r = ', error_r)
print('T = ', T, ',', 'M =', M, ',', 'Lambda =', Lambda, ',', 'K =', K, ',', 'niu_w = ', niu_w, ',', 'niu_f = ', niu_f, ',', 'miu = ', miu, ',', 'tau = ', tau, ',', 'beta = ', beta, ',', 'B = ', B, ',', 'k_f = ', k_f, ',', 'EPS = ', EPS)
# Polt
plt.figure()
plot(phi, title="phi")
phi_e_inp = interpolate(phi_e, P1Space)
plt.figure()
plot(phi_e_inp, title="phi_e")
plt.figure()
plot(p, title="p")
pe_inp = interpolate(p_e, P2Space)
plt.figure()
plot(pe_inp, title="pe")
plt.figure()
plot(u, title="u")
ue_inp = interpolate(u_e, P1_dgSpace)
plt.figure()
plot(ue_inp, title="ue")
plt.figure()
plot(m, title="m")
me_inp = interpolate(m_e, P1VectorSpace)
plt.figure()
plot(me_inp, title="me")
plt.figure()
plot(varphi, title="varphi")
varphie_inp = interpolate(varphi_e, P2Space)
plt.figure()
plot(varphie_inp, title="varphie")
plt.show()
In theory, the value of r and Sn should be similar, but after I run the program, the value of r and Sn is very different, and r is gradually decreasing. How to modify here ?
The following is the result when the running grid is divided into 1.0 / 128
r0= 1.4142135623730951
t= 0.0078125
r= 1.4142135623730951
S_n= 1.4142135623730951
t= 0.015625
r= 1.414187910672441
S_n= 1.4142058154844845
t= 0.0234375
r= 1.4141286188315734
S_n= 1.414179140332573
t= 0.03125
r= 1.414030554744989
S_n= 1.414131289019339
t= 0.0390625
r= 1.4138895806313179
S_n= 1.4140617500037005
t= 0.046875
r= 1.4137018768898324
S_n= 1.413970510188084
t= 0.0546875
r= 1.4134637320441845
S_n= 1.4138576606630635
t= 0.0625
r= 1.4131714905061583
S_n= 1.4137232960791615
t= 0.0703125
r= 1.4128215480407609
S_n= 1.4135674998836976
t= 0.078125
r= 1.4124103595575561
S_n= 1.4133903490259934
t= 0.0859375
r= 1.41193444907977
S_n= 1.413191920139061
t= 0.09375
r= 1.4113904195071811
S_n= 1.4129722937287472
t= 0.1015625
r= 1.410774961890038
S_n= 1.4127315565220093
t= 0.109375
r= 1.4100848643572623
S_n= 1.4124698026585278
t= 0.1171875
r= 1.4093170208441723
S_n= 1.4121871342497427
t= 0.125
r= 1.40846843969748
S_n= 1.4118836616154369
t= 0.1328125
r= 1.4075362521804875
S_n= 1.4115595033618282
t= 0.140625
r= 1.406517720867978
S_n= 1.4112147863828683
t= 0.1484375
r= 1.4054102479020651
S_n= 1.4108496458237658
t= 0.15625
r= 1.4042113830708816
S_n= 1.4104642250248425
t= 0.1640625
r= 1.4029188316673795
S_n= 1.4100586754540525
t= 0.171875
r= 1.4015304620834195
S_n= 1.4096331566318487
t= 0.1796875
r= 1.4000443130935962
S_n= 1.4091878360500254
t= 0.1875
r= 1.3984586007833693
S_n= 1.4087228890853771
t= 0.1953125
r= 1.396771725076741
S_n= 1.4082384989084076
t= 0.203125
r= 1.3949822758198587
S_n= 1.4077348563873093
t= 0.2109375
r= 1.3930890383784407
S_n= 1.407212159987287
t= 0.21875
r= 1.391090998708873
S_n= 1.4066706156653337
t= 0.2265625
r= 1.3889873478650303
S_n= 1.4061104367604265
t= 0.234375
r= 1.3867774859055253
S_n= 1.4055318438793059
t= 0.2421875
r= 1.384461025169007
S_n= 1.4049350647777654
t= 0.25
r= 1.382037792888353
S_n= 1.4043203342376596
t= 0.2578125
r= 1.379507833118138
S_n= 1.4036878939395057
t= 0.265625
r= 1.3768714079535125
S_n= 1.4030379923309333
t= 0.2734375
r= 1.3741289980226528
S_n= 1.4023708844909075
t= 0.28125
r= 1.3712813022391266
S_n= 1.401686831989864
t= 0.2890625
r= 1.3683292368049065
S_n= 1.400986102745791
t= 0.296875
r= 1.3652739334592163
S_n= 1.4002689708764173
t= 0.3046875
r= 1.362116736973009
S_n= 1.3995357165474651
t= 0.3125
r= 1.3588592018934875
S_n= 1.3987866258171575
t= 0.3203125
r= 1.355503088547701
S_n= 1.398021990477051
t= 0.328125
r= 1.3520503583188697
S_n= 1.3972421078892323
t= 0.3359375
r= 1.3485031682135893
S_n= 1.3964472808200579
t= 0.34375
r= 1.3448638647424989
S_n= 1.395637817270472
t= 0.3515625
r= 1.34113497714122
S_n= 1.3948140303030547
t= 0.359375
r= 1.337319209962443
S_n= 1.3939762378658813
t= 0.3671875
r= 1.3334194350738735
S_n= 1.3931247626133236
t= 0.375
r= 1.329438683100274
S_n= 1.3922599317239042
t= 0.3828125
r= 1.325380134351142
S_n= 1.391382076715329
t= 0.390625
r= 1.3212471092785079
S_n= 1.3904915332568084
t= 0.3984375
r= 1.317043058511857
S_n= 1.3895886409788183
t= 0.40625
r= 1.3127715525195578
S_n= 1.3886737432804
t= 0.4140625
r= 1.308436270947835
S_n= 1.3877471871341847
t= 0.421875
r= 1.3040409916899243
S_n= 1.3868093228892306
t= 0.4296875
r= 1.299589579739034
S_n= 1.3858605040718193
t= 0.4375
r= 1.2950859758793574
S_n= 1.384901087184382
t= 0.4453125
r= 1.290534185269625
S_n= 1.3839314315027063
t= 0.453125
r= 1.2859382659735514
S_n= 1.3829518988715321
t= 0.4609375
r= 1.2813023174908689
S_n= 1.3819628534987236
t= 0.46875
r= 1.2766304693418318
S_n= 1.3809646617481495
t= 0.4765625
r= 1.2719268697566168
S_n= 1.3799576919314775
t= 0.484375
r= 1.2671956745195159
S_n= 1.378942314098959
t= 0.4921875
r= 1.2624410360157536
S_n= 1.377918899829465
t= 0.5
r= 1.2576670925265407
S_n= 1.3768878220198562