我的 durand-kerner 实现有什么问题?
What's wrong with my durand-kerner implementation?
实现这个简单的寻根算法。http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method我一辈子都无法弄清楚我的实施出了什么问题。根部不断爆炸,没有收敛的迹象。有什么建议吗?
谢谢。
#include <iostream>
#include <complex>
using namespace std;
typedef complex<double> dcmplx;
dcmplx f(dcmplx x)
{
// the function we are interested in
double a4 = 3;
double a3 = -3;
double a2 = 1;
double a1 = 0;
double a0 = 100;
return a4 * pow(x,4) + a3 * pow(x,3) + a2 * pow(x,2) + a1 * x + a0;
}
int main()
{
dcmplx p(.9,2);
dcmplx q(.1, .5);
dcmplx r(.7,1);
dcmplx s(.3, .5);
dcmplx p0, q0, r0, s0;
int max_iterations = 20;
bool done = false;
int i=0;
while (i<max_iterations && done == false)
{
p0 = p;
q0 = q;
r0 = r;
s0 = s;
p = p0 - f(p0)/((p0-q0)*(p0-r0)*(p0-s0));
q = q0 - f(q0)/((q0-p)*(q0-r0)*(q0-s0));
r = r0 - f(r0)/((r0-p)*(r0-q)*(r0-s0));
s = s0 - f(s0)/((s0-p)*(s0-q)*(s0-r));
// if convergence within small epsilon, declare done
if (abs(p-p0)<1e-5 && abs(q-q0)<1e-5 && abs(r-r0)<1e-5 && abs(s-s0)<1e-5)
done = true;
i++;
}
cout<<"roots are :n";
cout << p << "n";
cout << q << "n";
cout << r << "n";
cout << s << "n";
cout << "number steps taken: "<< i << endl;
return 0;
}
晚了半年:谜团的解决方案是分母应该是多项式导数的近似值,因此需要包含前导系数 a4 作为因子。
或者,可以在 return 语句中将多项式值除以 a4,以便多项式有效赋范,即具有前导系数 1。
请注意,Bo Jacoby 在维基百科中的示例代码是该方法的 Seidel 型变体,经典公式是类似乔丹的方法,其中所有新的近似值都是从旧的近似同时计算出来的。Seidel可以具有比作为多维牛顿方法的公式为雅可比提供的阶2更快的收敛。
然而,对于大度雅可比,可以使用快速多项式乘法算法加速,以对多项式值和分母中的乘积进行所需的多点评估。
啊,问题是 N 次多项式的系数必须指定为
1*x^N + a*x^(N-1) + b*x^(N-2) ...等 + z;
其中 1 是最大度项的系数。否则,第一个根将永远不会收敛。
您尚未正确实现公式。例如
s = s0 - f(s0)/((s0-p0)*(s0-q0)*(s0-r0));
应该是
s = s0 - f(s0)/((s0-p)*(s0-q)*(s0-r));
再看维基文章
相关文章:
- 警告处理为错误这里有什么问题
- C++我的数学有什么问题,为什么我的代码不能正确循环
- 问题:什么是QAbstractItemView::NoEditTriggers的反面
- 当我尝试添加 2 个大字符串时,我无法弄清楚出了什么问题
- 违反const正确性:我应该现实地期待什么问题
- 这个带有模板<类 Vector 的C++代码片段有什么问题>
- 我的逻辑反转字符串中的元音有什么问题?
- 需要以下代码的帮助,下面的代码有什么问题
- 常量公共成员有什么问题?
- c++无值sort()的问题是什么?
- 以下代码中的函数模板有什么问题?
- 这个返回元素位置的基于循环的函数有什么问题?
- creat_list2功能有什么问题?
- 基本的 c++ 问题:如果我在函数中创建某些内容并返回它会发生什么?
- 我遇到了黑客排名中的问题"TWO STRINGS"的三个测试用例的分段错误。原因是什么?
- 什么是钻石问题?是一系列问题还是特定问题?
- 格式说明符C++有什么问题
- 我应该在 main 函数中写什么来测试我的问题?
- 任何人都可以告诉我我的 C++ 代码出了什么问题?
- 方法问题 - 什么会改变值,什么不会改变?什么是无效的?