求勒让德多项式的递归程序
Recursive program to obtain Legendre Polynomials
我正在尝试在 C++ 上为 Legendre 多项式的构造函数编写一个函数,该函数打印出次数为 m 的多项式的系数。多项式遵循 simple recursion relation.
现在,我正在尝试实现该关系,并且它对每个 n 都完美工作,最多 8 个,从 9 开始,在第 7 次迭代周围,它会拾取一个垃圾值甚至不在任何先前的系数向量中。我想知道我该如何解决这个问题。我向您展示我的代码:
#include <vector>
#include <cmath>
std::vector<double> set_coeffs(int m){
std::vector<double> coeffs;
if (m == 0) //Casos iniciales para empezar la recursión
{
coeffs.push_back(1);
} else if (m == 1)
{
coeffs.push_back(1);
coeffs.push_back(0);
} else if (m == 2) //Puse también el caso 2 porque de otro modo ocurre el mismo problema pero desde n=5
{
coeffs.push_back(1.5);
coeffs.push_back(0);
coeffs.push_back(-0.5);
}
else
{
std::vector<double> v = set_coeffs(m-1);
std::vector<double> u = set_coeffs(m-2);
std::cout << "inicia cicle" << std::endl;
double a = (2* ((double)m) -1)/((double)m);
double b = (((double)m)-1)/((double)m);
coeffs.push_back(a*v[0]);
coeffs.push_back(a*v[1]);
for (int i = 0; i < m-1; i++)
{
double c = a*v[i+2] - b*u[i];
std::cout << m << " " << v[i+2] << " " << u[i] << " " << c <<std::endl;
coeffs.push_back(c);
}
std::cout << "termina cicle" << std::endl;
}
return coeffs;
}
正在调用 set_coeffs(n)
returns 具有 n+1 个元素的向量。
因此,std::vector<double> v = set_coeffs(m-1);
有 m 个元素。
在您的 for
循环中,i
从 0
运行到 m-2
,然后您访问 v[i+2]
。在最后一次迭代中,这将访问超出范围的 v[m]
。
如 interjay's 中所述,发布的代码没有正确考虑向量大小的差异。
为了使索引匹配更容易(至少恕我直言,但它也应该有助于涉及这些系数的进一步计算),我将以相反的顺序存储系数。我也会利用它们的特定模式(只有奇数 或 偶数不为零)。
可能值得一提的是,当应用于大向量时,这种递归方法效率低下。
也就是说,这是一个可能的实现:
std::vector<double> Legendre_coefficients(int m)
{
if (m == 0)
{
return {1};
}
if (m == 1)
{
return {0, 1};
}
// Initialize with zero, only at most (half + 1) of the terms will be changed later
std::vector<double> coeffs(m + 1);
// Consider some form of memoization instead of this recursion
std::vector<double> v = Legendre_coefficients(m - 1);
std::vector<double> u = Legendre_coefficients(m - 2);
// using literals of floating point type, 'm' is promoted by the compiler
double a = (2.0 * m - 1.0) / m;
double b = (m - 1.0) / m;
int first = 1;
// If 'm' is even, so is (m - 2) and its first element is zero. It can be skipped.
// It also avoids annoying '-0' in the output
if ( m % 2 == 0 )
{
coeffs[0] = -b * u[0];
first = 2;
}
for (int i = first; i < m - 1; i += 2)
{
coeffs[i] = (a * v[i - 1] - b * u[i]);
}
coeffs[m] = a * v[m - 1];
return coeffs;
}
可测试here。
我正在尝试在 C++ 上为 Legendre 多项式的构造函数编写一个函数,该函数打印出次数为 m 的多项式的系数。多项式遵循 simple recursion relation.
现在,我正在尝试实现该关系,并且它对每个 n 都完美工作,最多 8 个,从 9 开始,在第 7 次迭代周围,它会拾取一个垃圾值甚至不在任何先前的系数向量中。我想知道我该如何解决这个问题。我向您展示我的代码:
#include <vector>
#include <cmath>
std::vector<double> set_coeffs(int m){
std::vector<double> coeffs;
if (m == 0) //Casos iniciales para empezar la recursión
{
coeffs.push_back(1);
} else if (m == 1)
{
coeffs.push_back(1);
coeffs.push_back(0);
} else if (m == 2) //Puse también el caso 2 porque de otro modo ocurre el mismo problema pero desde n=5
{
coeffs.push_back(1.5);
coeffs.push_back(0);
coeffs.push_back(-0.5);
}
else
{
std::vector<double> v = set_coeffs(m-1);
std::vector<double> u = set_coeffs(m-2);
std::cout << "inicia cicle" << std::endl;
double a = (2* ((double)m) -1)/((double)m);
double b = (((double)m)-1)/((double)m);
coeffs.push_back(a*v[0]);
coeffs.push_back(a*v[1]);
for (int i = 0; i < m-1; i++)
{
double c = a*v[i+2] - b*u[i];
std::cout << m << " " << v[i+2] << " " << u[i] << " " << c <<std::endl;
coeffs.push_back(c);
}
std::cout << "termina cicle" << std::endl;
}
return coeffs;
}
正在调用 set_coeffs(n)
returns 具有 n+1 个元素的向量。
因此,std::vector<double> v = set_coeffs(m-1);
有 m 个元素。
在您的 for
循环中,i
从 0
运行到 m-2
,然后您访问 v[i+2]
。在最后一次迭代中,这将访问超出范围的 v[m]
。
如 interjay's
为了使索引匹配更容易(至少恕我直言,但它也应该有助于涉及这些系数的进一步计算),我将以相反的顺序存储系数。我也会利用它们的特定模式(只有奇数 或 偶数不为零)。
可能值得一提的是,当应用于大向量时,这种递归方法效率低下。
也就是说,这是一个可能的实现:
std::vector<double> Legendre_coefficients(int m)
{
if (m == 0)
{
return {1};
}
if (m == 1)
{
return {0, 1};
}
// Initialize with zero, only at most (half + 1) of the terms will be changed later
std::vector<double> coeffs(m + 1);
// Consider some form of memoization instead of this recursion
std::vector<double> v = Legendre_coefficients(m - 1);
std::vector<double> u = Legendre_coefficients(m - 2);
// using literals of floating point type, 'm' is promoted by the compiler
double a = (2.0 * m - 1.0) / m;
double b = (m - 1.0) / m;
int first = 1;
// If 'm' is even, so is (m - 2) and its first element is zero. It can be skipped.
// It also avoids annoying '-0' in the output
if ( m % 2 == 0 )
{
coeffs[0] = -b * u[0];
first = 2;
}
for (int i = first; i < m - 1; i += 2)
{
coeffs[i] = (a * v[i - 1] - b * u[i]);
}
coeffs[m] = a * v[m - 1];
return coeffs;
}
可测试here。