From 75668b8e561c301c9d1fc9cd37e2dc487dd1c2d4 Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Thu, 7 Mar 2024 10:37:58 +0000 Subject: [PATCH] [ode] Remove unused test (amplifier) --- ode/problems.go | 104 +++------------------------------------------- ode/t_ode_test.go | 40 ------------------ 2 files changed, 6 insertions(+), 138 deletions(-) diff --git a/ode/problems.go b/ode/problems.go index 21394a2c..5ff37e0c 100644 --- a/ode/problems.go +++ b/ode/problems.go @@ -65,7 +65,8 @@ func (o *Problem) Solve(method string, fixedStp, numJac bool) (y la.Vector, stat } // ConvergenceTest runs convergence test -// yExact -- is the exact (reference) y @ xf +// +// yExact -- is the exact (reference) y @ xf func (o *Problem) ConvergenceTest(tst *testing.T, dxmin, dxmax float64, ndx int, yExact la.Vector, methods []string, orders, tols []float64) { @@ -145,9 +146,10 @@ func ProbHwEq11() (o *Problem) { } // ProbVanDerPol returns the Van der Pol' Equation as given in Hairer-Wanner VII-p5 Eq.(1.5) -// eps -- ε coefficient; use 0 for default value [=1e-6] -// stationary -- use eps=1 and compute period and amplitude such that -// y = [A, 0] is a stationary point +// +// eps -- ε coefficient; use 0 for default value [=1e-6] +// stationary -- use eps=1 and compute period and amplitude such that +// y = [A, 0] is a stationary point func ProbVanDerPol(eps float64, stationary bool) (o *Problem) { o = new(Problem) @@ -218,100 +220,6 @@ func ProbRobertson() (o *Problem) { return } -// ProbHwAmplifier returns the Hairer-Wanner VII-p376 Transistor Amplifier -// NOTE: from E. Hairer's website, not the equation in the book -func ProbHwAmplifier() (o *Problem) { - - // problem - o = new(Problem) - o.Xf = 0.05 - //o.Xf = 0.0123 // OK - //o.Xf = 0.01235 // !OK - - // constants - ue, ub, uf, α, β := 0.1, 6.0, 0.026, 0.99, 1.0e-6 - r0, r1, r2, r3, r4, r5 := 1000.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0 - r6, r7, r8, r9 := 9000.0, 9000.0, 9000.0, 9000.0 - w := 2.0 * 3.141592654 * 100.0 - - // initial values - o.Y = la.NewVectorSlice([]float64{ - 0.0, - ub, - ub / (r6/r5 + 1.0), - ub / (r6/r5 + 1.0), - ub, - ub / (r2/r1 + 1.0), - ub / (r2/r1 + 1.0), - 0.0, - }) - o.Ndim = len(o.Y) - - // right-hand side of the amplifier problem - o.Fcn = func(f la.Vector, dx, x float64, y la.Vector) { - uet := ue * math.Sin(w*x) - fac1 := β * (math.Exp((y[3]-y[2])/uf) - 1.0) - fac2 := β * (math.Exp((y[6]-y[5])/uf) - 1.0) - f[0] = y[0] / r9 - f[1] = (y[1]-ub)/r8 + α*fac1 - f[2] = y[2]/r7 - fac1 - f[3] = y[3]/r5 + (y[3]-ub)/r6 + (1.0-α)*fac1 - f[4] = (y[4]-ub)/r4 + α*fac2 - f[5] = y[5]/r3 - fac2 - f[6] = y[6]/r1 + (y[6]-ub)/r2 + (1.0-α)*fac2 - f[7] = (y[7] - uet) / r0 - } - - // Jacobian of the amplifier problem - o.Jac = func(dfdy *la.Triplet, dx, x float64, y la.Vector) { - fac14 := β * math.Exp((y[3]-y[2])/uf) / uf - fac27 := β * math.Exp((y[6]-y[5])/uf) / uf - if dfdy.Max() == 0 { - dfdy.Init(8, 8, 16) - } - nu := 2 - dfdy.Start() - dfdy.Put(2+0-nu, 0, 1.0/r9) - dfdy.Put(2+1-nu, 1, 1.0/r8) - dfdy.Put(1+2-nu, 2, -α*fac14) - dfdy.Put(0+3-nu, 3, α*fac14) - dfdy.Put(2+2-nu, 2, 1.0/r7+fac14) - dfdy.Put(1+3-nu, 3, -fac14) - dfdy.Put(2+3-nu, 3, 1.0/r5+1.0/r6+(1.0-α)*fac14) - dfdy.Put(3+2-nu, 2, -(1.0-α)*fac14) - dfdy.Put(2+4-nu, 4, 1.0/r4) - dfdy.Put(1+5-nu, 5, -α*fac27) - dfdy.Put(0+6-nu, 6, α*fac27) - dfdy.Put(2+5-nu, 5, 1.0/r3+fac27) - dfdy.Put(1+6-nu, 6, -fac27) - dfdy.Put(2+6-nu, 6, 1.0/r1+1.0/r2+(1.0-α)*fac27) - dfdy.Put(3+5-nu, 5, -(1.0-α)*fac27) - dfdy.Put(2+7-nu, 7, 1.0/r0) - } - - // "mass" matrix - c1, c2, c3, c4, c5 := 1.0e-6, 2.0e-6, 3.0e-6, 4.0e-6, 5.0e-6 - o.M = new(la.Triplet) - o.M.Init(8, 8, 14) - o.M.Start() - nu := 1 - o.M.Put(1+0-nu, 0, -c5) - o.M.Put(0+1-nu, 1, c5) - o.M.Put(2+0-nu, 0, c5) - o.M.Put(1+1-nu, 1, -c5) - o.M.Put(1+2-nu, 2, -c4) - o.M.Put(1+3-nu, 3, -c3) - o.M.Put(0+4-nu, 4, c3) - o.M.Put(2+3-nu, 3, c3) - o.M.Put(1+4-nu, 4, -c3) - o.M.Put(1+5-nu, 5, -c2) - o.M.Put(1+6-nu, 6, -c1) - o.M.Put(0+7-nu, 7, c1) - o.M.Put(2+6-nu, 6, c1) - o.M.Put(1+7-nu, 7, -c1) - return -} - // ProbArenstorf returns the Arenstorf orbit problem func ProbArenstorf() (o *Problem) { o = new(Problem) diff --git a/ode/t_ode_test.go b/ode/t_ode_test.go index 66c5de94..116c2fda 100644 --- a/ode/t_ode_test.go +++ b/ode/t_ode_test.go @@ -6,7 +6,6 @@ package ode import ( "testing" - "time" "github.com/cpmech/gosl/chk" "github.com/cpmech/gosl/io" @@ -149,42 +148,3 @@ func TestOde03(tst *testing.T) { chk.Int(tst, "number of lin solutions ", sol.Stat.Nlinsol, 24) chk.Int(tst, "max number of iterations", sol.Stat.Nitmax, 2) } - -func TestOde04(tst *testing.T) { - - //verbose() - chk.PrintTitle("ode04: Hairer-Wanner VII-p376 Transistor Amplifier") - - // problem - p := ProbHwAmplifier() - - // configurations - conf := NewConfig("radau5", "") - conf.SetStepOut(true, nil) - conf.IniH = 1.0e-6 // initial step size - - // set tolerances - atol, rtol := 1e-11, 1e-5 - conf.SetTols(atol, rtol) - - // ODE solver - sol := NewSolver(p.Ndim, conf, p.Fcn, p.Jac, p.M) - defer sol.Free() - - // run - t0 := time.Now() - sol.Solve(p.Y, 0.0, p.Xf) - io.Pfmag("elapsed time = %v\n", time.Now().Sub(t0)) - - // check - if false { // these values vary slightly in different machines - chk.Int(tst, "number of F evaluations ", sol.Stat.Nfeval, 2599) - chk.Int(tst, "number of J evaluations ", sol.Stat.Njeval, 216) - chk.Int(tst, "total number of steps ", sol.Stat.Nsteps, 275) - chk.Int(tst, "number of accepted steps", sol.Stat.Naccepted, 219) - chk.Int(tst, "number of rejected steps", sol.Stat.Nrejected, 20) - chk.Int(tst, "number of decompositions", sol.Stat.Ndecomp, 274) - chk.Int(tst, "number of lin solutions ", sol.Stat.Nlinsol, 792) - chk.Int(tst, "max number of iterations", sol.Stat.Nitmax, 6) - } -}