Golang Environment.Beta方法代码示例

// Partial derivative of F with respect to T; Mu_h and V held constant.
func dFdT(env *tempAll.Environment, F envFunc) (float64, error) {
	ct := 0
	// G gets F given Beta (allow x to vary; constant Mu_h)
	G := func(Beta float64) (float64, error) {
		ct += 1
		// save the environment state before changing it
		// (don't want one call of F to affect the next)
		oD1, oBeta, oX, oMu_b := env.D1, env.Beta, env.X, env.Mu_b
		env.Beta = Beta
		// fix free variables
		eps := 1e-9
		_, err := D1F0XSolve(env, eps, eps)
		if err != nil {
			return 0.0, err
		}
		vF, err := F(env)
		if err != nil {
			return 0.0, err
		}
		// restore the environment
		env.D1, env.Beta, env.X, env.Mu_b = oD1, oBeta, oX, oMu_b
		return vF, nil
	}
	h := 1e-4
	epsAbs := 1e-5
	deriv, err := solve.OneDimDerivative(G, env.Beta, h, epsAbs)
	//fmt.Println("dF_dT ct", ct)
	return -math.Pow(env.Beta, 2.0) * deriv, err
}

// Partial derivative of Mu_h with respect to T; x and V held constant.
func dMu_hdT(env *tempAll.Environment) (float64, error) {
	ct := 0
	// F gets Mu_h given Beta
	F := func(Beta float64) (float64, error) {
		ct += 1
		// save the environment state before changing it
		// (don't want one call of F to affect the next)
		oD1, oMu_h, oBeta, oMu_b := env.D1, env.Mu_h, env.Beta, env.Mu_b
		env.Beta = Beta
		// fix free variables
		eps := 1e-9
		_, err := D1MuF0Solve(env, eps, eps)
		if err != nil {
			return 0.0, err
		}
		Mu_h := env.Mu_h
		// restore the environment
		env.D1, env.Mu_h, env.Beta, env.Mu_b = oD1, oMu_h, oBeta, oMu_b
		return Mu_h, nil
	}
	h := 1e-4
	epsAbs := 1e-5
	deriv, err := solve.OneDimDerivative(F, env.Beta, h, epsAbs)
	//fmt.Println("MuT ct", ct)
	return -math.Pow(env.Beta, 2.0) * deriv, err
}

// Solve the (D1, Mu_h, Beta) system with x and F0 fixed.
func D1MuBetaSolve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	// our guess for beta should be above beta_c
	if env.A == 0.0 && env.B == 0.0 {
		D1, Mu_h, F0 := env.D1, env.Mu_h, env.F0
		env.F0 = 0.0 // F0 is 0 at T_c
		_, err := tempCrit.CritTempSolve(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		fmt.Printf("%v; Tc = %f\n", env, 1.0/env.Beta)
		omegaFit, err := tempCrit.OmegaFit(env, tempCrit.OmegaPlus)
		if err != nil {
			return nil, err
		}
		env.A, env.B = omegaFit[0], omegaFit[2]
		env.PairCoeffsReady = true
		env.Beta += 0.1
		// we are at T < T_c; uncache env
		env.D1, env.Mu_h, env.F0 = D1, Mu_h, F0
	}
	//fmt.Printf("%v; Tc = %f\n", env, 1.0 / env.Beta)
	// solve low temp system for reasonable values of D1 and Mu_h first
	_, err := D1MuSolve(env, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full low temp system
	system, start := D1MuBetaSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

func SolveNoninteracting(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	env.F0 = 0.0
	env.Mu_h = 0.3
	env.Beta = 50.0
	system, start := NoninteractingSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

// For use with solve.Iterative:
func CritTempStages(env *tempAll.Environment) ([]solve.DiffSystem, []vec.Vector, func([]vec.Vector)) {
	vars0 := []string{"D1", "Mu_h"}
	vars1 := []string{"Beta"}
	diffD1 := tempPair.AbsErrorD1(env, vars0)
	diffMu_h := tempPair.AbsErrorBeta(env, vars0)
	system0 := solve.Combine([]solve.Diffable{diffD1, diffMu_h})
	diffBeta := AbsErrorBeta(env, vars1)
	system1 := solve.Combine([]solve.Diffable{diffBeta})
	stages := []solve.DiffSystem{system0, system1}
	start := []vec.Vector{[]float64{env.D1, env.Mu_h}, []float64{env.Beta}}
	accept := func(x []vec.Vector) {
		env.D1 = x[0][0]
		env.Mu_h = x[0][1]
		env.Beta = x[1][0]
	}
	return stages, start, accept
}

// Solve the environment under the conditions at T = T_c.
func CritTempSolve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	// our guess for beta should be a bit above Beta_p
	pairSystem, pairStart := tempPair.PairTempSystem(env)
	_, err := solve.MultiDim(pairSystem, pairStart, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	env.Beta += 0.1
	// solve crit temp system for reasonable values of Mu and D1 first
	system, start := CritTempD1MuSystem(env)
	_, err = solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full crit temp system
	system, start = CritTempFullSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}

// Solve the (D1, Mu_h, Beta) system with x and Mu_b fixed.
func FlucTempSolve(env *tempAll.Environment, epsAbs, epsRel float64) (vec.Vector, error) {
	// fix pair coefficients
	if env.A == 0.0 && env.B == 0.0 && env.FixedPairCoeffs {
		D1, Mu_h, Mu_b, Beta := env.D1, env.Mu_h, env.Mu_b, env.Beta
		env.Mu_b = 0.0 // Mu_b is 0 at T_c
		_, err := tempCrit.CritTempSolve(env, epsAbs, epsRel)
		if err != nil {
			return nil, err
		}
		omegaFit, err := tempCrit.OmegaFit(env, tempCrit.OmegaPlus)
		if err != nil {
			return nil, err
		}
		env.A, env.B = omegaFit[0], omegaFit[2]
		env.PairCoeffsReady = true
		// uncache env
		env.D1, env.Mu_h, env.Mu_b, env.Beta = D1, Mu_h, Mu_b, Beta
	}
	// our guess for beta should be a bit above Beta_p
	pairSystem, pairStart := tempPair.PairTempSystem(env)
	_, err := solve.MultiDim(pairSystem, pairStart, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	env.Beta += 0.1
	// solve fluc temp system for reasonable values of Mu_h and D1 first
	system, start := FlucTempD1MuSystem(env)
	_, err = solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	// solve the full fluc temp system
	system, start = FlucTempFullSystem(env)
	solution, err := solve.MultiDim(system, start, epsAbs, epsRel)
	if err != nil {
		return nil, err
	}
	return solution, nil
}