Skip to content

Commit

Permalink
[ode] Remove unused test (amplifier)
Browse files Browse the repository at this point in the history
  • Loading branch information
cpmech committed Mar 7, 2024
1 parent a62d39e commit 75668b8
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 138 deletions.
104 changes: 6 additions & 98 deletions ode/problems.go
Original file line number Diff line number Diff line change
Expand Up @@ -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) {

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
40 changes: 0 additions & 40 deletions ode/t_ode_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ package ode

import (
"testing"
"time"

"github.com/cpmech/gosl/chk"
"github.com/cpmech/gosl/io"
Expand Down Expand Up @@ -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)
}
}

0 comments on commit 75668b8

Please sign in to comment.