variational Problem - reviewed changes

Author: arnavk23

src/ADNLPProblems/variational.jl

@@ -2,25 +2,30 @@ export variational
 
 function variational(; n::Int = default_nvar, type::Type{T} = Float64, kwargs...) where {T}
   h = 1 // (n + 1)
-  x0 = Vector{T}(undef, n)
-  for i = 1:n
-    x0[i] = i * h * (1 - i * h)
-  end
-  f =
-    x -> begin
-      xext = vcat(zero(eltype(x)), x, zero(eltype(x)))
-      term1 = sum(x[i] * (x[i] - xext[i + 2]) / convert(eltype(x), h) for i = 1:n)
-      term2 = sum(
-        if eltype(x) <: AbstractFloat
-          d = xext[j + 2] - xext[j + 1]
-          abs(d) <= eps(eltype(x)) ? exp(xext[j + 1]) : (exp(xext[j + 2]) - exp(xext[j + 1])) / d
-        else
-          d = xext[j + 2] - xext[j + 1]
-          exp(xext[j + 1]) * expm1(d) / d
-        end for j = 0:n
-      )
-      return 2 * (term1 + n * (convert(eltype(x), h) / 2) * term2)
+  x0 = [
+    begin
+      ih = i * h
+      convert(T, ih * (1 - ih))
+    end for i = 1:n
+  ]
+
+  function f(x)
+    term1 = zero(T)
+    for i = 1:n
+      xi = x[i]
+      xip = (i < n) ? x[i + 1] : zero(T)
+      term1 += xi * (xi - xip) / h
+    end
+
+    term2 = zero(T)
+    for j = 0:n
+      a = (j == 0) ? zero(T) : x[j]
+      b = (j == n) ? zero(T) : x[j + 1]
+      term2 += (exp(b) - exp(a)) / (b - a)
     end
 
+    return 2 * (term1 + n * (h / 2) * term2)
+  end
+
   return ADNLPModels.ADNLPModel(f, x0, name = "variational"; kwargs...)
 end

src/PureJuMP/variational.jl

@@ -12,10 +12,12 @@ export variational
 
 function variational(; n::Int = default_nvar, kwargs...)
   h = 1 // (n + 1)
-    x0 = [begin
-        ih = i * h
-        convert(Float64, ih * (1 - ih))
-      end for i = 1:n]
+  x0 = [
+    begin
+      ih = i * h
+      convert(T, ih * (1 - ih))
+    end for i = 1:n
+  ]
   model = Model()
   @variable(model, x[i = 1:n], start = x0[i])
 
@@ -27,17 +29,9 @@ function variational(; n::Int = default_nvar, kwargs...)
       n *
       (h / 2) *
       (
-        (1 + x[1]/2 + x[1]^2/6 + x[1]^3/24 + x[1]^4/120) +
-        sum(
-          exp(x[j]) * (
-            1 +
-            (x[j + 1] - x[j]) / 2 +
-            (x[j + 1] - x[j])^2 / 6 +
-            (x[j + 1] - x[j])^3 / 24 +
-            (x[j + 1] - x[j])^4 / 120
-          ) for j = 1:(n - 1)
-        ) +
-        exp(x[n]) * (1 + (0 - x[n]) / 2 + (0 - x[n])^2 / 6 + (0 - x[n])^3 / 24 + (0 - x[n])^4 / 120)
+        (exp(x[1]) - exp(0)) / (x[1] - 0) +
+        sum((exp(x[j + 1]) - exp(x[j])) / (x[j + 1] - x[j]) for j = 1:(n - 1)) +
+        (exp(0) - exp(x[n])) / (0 - x[n])
       )
     )
   )