1

假设我们有一个二维函数 f(x,y)，我们想用一组切比雪夫多项式进行逼近，直到 2 次。设 j 次切比雪夫多项式为 T j (x) 或 T j (y)。我们通常通过构造一个函数 g(x,y) 来近似 f(x,y)，它是一维多项式的张量积，

        

我想要做的是生成一个完整的 N 级切比雪夫多项式。这只是上面的张量积，但是索引 k+l 的总和必须小于或等于 N。所以如果 N 为 3，那么我们将得到除 T 2 (x)*T 2 (y) 之外的所有项，因为 2+2=4 > 3。随着函数维数的增加，更多项会被删除。

最终，如果我使用超过 2 或 3 个维度，我希望通过灵活的级别选择来做到这一点，并且不必写出一堆嵌套循环。这似乎@nloops是要走的路，但我无法弄清楚。

例如，假设我想在 (.5,.5) 处评估二维切比雪夫多项式。我可以编写一个内联函数，它在点 x 处返回 N 级的一维切比雪夫多项式。

cheb(d,x) = cos((d)*acos(x))
a = [.5, .5]
polys1d = [cheb(d,a[i]) for d = 0:2, i = 1:length(a)]

在 2 维（甚至更多）中创建完整的张量积多项式很容易。例如：

polys2d = kron(polys1d[:,1],polys1d[:,2])

但是以一般方式创建完整的多项式有点棘手。我希望它以一种我不通过构造完整张量积然后删除度数总和大于水平的多项式开始的方式来完成。如果维度的数量和级别都很大，这将占用大量的内存。

4

1 回答 1

1

如果我们定义以下辅助函数：

using Iterators   # install with Pkg.add("Iterators")

nlimitedkparts(n,k) = (diff([0;v]).-1 for v in subsets(1:(n+k),k))

我们可以生成以下内容：

julia> ["T_$(r[1])(x)*T_$(r[2])(y)" for r in nlimitedkparts(3,2)]
10-element Array{String,1}:
 "T_0(x)*T_0(y)"
 "T_0(x)*T_1(y)"
 "T_0(x)*T_2(y)"
 "T_0(x)*T_3(y)"
 "T_1(x)*T_0(y)"
 "T_1(x)*T_1(y)"
 "T_1(x)*T_2(y)"
 "T_2(x)*T_0(y)"
 "T_2(x)*T_1(y)"
 "T_3(x)*T_0(y)"

当然，除了生成这些字符串之外，您可能还想做其他事情，但使用相同的nlimitedkparts函数。

于 2017-07-27T05:40:35.947 回答