0

我已经编写了一种通过复合辛普森规则逼近定积分的方法。

#=
f  integrand
a  lower integration bound
b  upper integration bound
n  number of iterations or panels

h     step size
=#
function simpson(f::Function, a::Number, b::Number, n::Number)
   n % 2 == 0 || error("`n` must be even")
   h = (b - a) / n
   s = f(a) + f(b)
   s += 4*sum(f(a .+ collect(1:2:n) .* h))
   s += 2*sum(f(a .+ collect(2:2:n-1) .* h))
   return h/3 * s
end

对于“简单”功能，例如e^(-x^2)，该simpson功能有效。

Input: simpson(x -> simpson(x -> exp.(-x.^2), 0, 5, 100)
Output: 0.8862269254513949

但是，对于更复杂的功能f(x)

gArgs(x) = (30 .+ x, 0)
f(x) = exp.(-x.^2) .* maximum(generator.(gArgs.(x)...)[1])

其中generator(θ, plotsol)是一个函数，它以百分比表示的缺陷 θ 和一个布尔值plotsol（0 或 1）确定是否应绘制发电机，并返回一个矢量，其中包含发电机中某些点的磁化强度。

当我尝试通过运行以下代码来计算积分时

gArgs(x) = (30 .+ x, 0)
f(x) = exp.(-x.^2) .* maximum(generator.(gArgs.(x)...)[1])
println(simpson(x -> f(x), 0, 5, 10))

我遇到错误MethodError: no method matching generator(::Float64)。对于表达式的轻微变体，f(x)我会遇到不同的错误，例如DimensionMismatch("array could not be broadcast to match destination")和InexactError: Bool(33.75)。最后，我认为错误的原因归结为我无法弄清楚如何正确输入被积函数的表达式f(x)。有人可以帮我弄清楚如何f(x)正确输入吗？如果我的问题有什么不清楚的地方，请告诉我。

4

1 回答 1

1

给定一个数组x，gArgs.(x)返回一个Tuples 数组，并且您正试图通过一个元组数组进行广播。但是使用元组广播的行为有点不同。元组不被视为单个元素，它们本身是广播的。

julia> println.(gArgs.([0.5, 1.5, 2.5, 3.5, 4.5])...)
30.531.532.533.534.5
00000

这不是你所期望的，是吗？

您还可以通过以下示例看到问题；

julia> (2, 5) .!= [(2, 5)]
2-element BitArray{1}:
 true
 true

我相信f这是一个实际上需要一个标量并返回一个标量的函数。f您应该将广播留给调用者，而不是处理数组。您很可能会更好地实现f元素。这是更加 Julia 的做事方式，会让你的工作更轻松。

也就是说，我相信您的实现应该与以下修改一起使用，如果您在 generator.

function simpson(f::Function, a::Number, b::Number, n::Number)
   n % 2 == 0 || error("`n` must be even")
   h = (b - a) / n
   s = f(a) + f(b)
   s += 4*sum(f.(a .+ collect(1:2:n) .* h)) # broadcast `f`
   s += 2*sum(f.(a .+ collect(2:2:n-1) .* h)) # broadcast `f`
   return h/3 * s
end

# define `gArg` and `f` element-wise and `generator`, too.
gArgs(x) = (30 + x, 0) # get rid of broadcasting dot. Shouldn't `0` be `false`?
f(x) = exp(-x^2) * maximum(generator(gArgs(x)...)[1]) # get rid of broadcasting dots
println(simpson(f, 0, 5, 10)) # you can just write `f`

您还应该按generator元素定义函数。

于 2019-04-17T16:49:29.270 回答