h**********c
发帖数: 4120 | 1 先前提到过 y = sin(a * pi /x), 对于质数a 在（1，a）是没有解得。
写了段程序，用牛顿法算 x， 跑了跑发现，原来这个一维系统有很多folds，又想了想
能不能用二维的系统来解这个问题呢？于是又了如下的构造
f = (f1,f2)^t
f1(x1,x2) = (x1 + exp(x2))* sin(pi * a /x1);
f2(x1,x2) = (x1*x1 -x2*x2)
求 f --〉0
这个系统看起来有个不错的jacobian，跑了跑程序，居然收敛了，但不知道收敛到哪里
去了，您有时间帮忙看看吧，能不能有更好的构造。或者用什么论证着干脆就是不可行
的，免得骑自行车去月球。
附程序c++
#include
#include
using namespace std;
#define PI (atan(1.0)*4.0)
#define __VERYSMALL (1.0E-10)
#define __ITERATIONLIMIT 20
bool isVerySmall (double d) {
return (abs(d)<__VERYSMALL);
}
// exp(x1+x2) is too big
inline double f1 (double x1, double x2,double a) {
//return exp(x1+x2) * sin( PI * a / x1);
return (x1+exp(x2)) * sin( PI * a / x1);
}
inline double f2 (double x1, double x2,double a) {
return x1*x1 - x2* x2;
}
inline double f1px1 (double x1, double x2,double a) {
return x1 * sin( PI * a / x1) - (x1 + exp(x2) )* cos(PI * a / x1) * PI *
a * (1.0/x1/x1) * exp(x1+x2) ;
}
inline double f1px2 (double x1, double x2,double a) {
//return exp(x1+x2) * sin( PI * a / x1);
return exp(x2) * sin( PI * a / x1);
}
inline double f2px1 (double x1, double x2,double a) {
return 2.0*x1;
}
inline double f2px2 (double x1, double x2,double a) {
return -2.0*x2;
}
void computeDeltaX (double * arr, double * f, double * result) {
double det = (arr[1]*arr[2] -arr[0]*arr[3]);
if (isVerySmall (det)) {
cout<< "det is very samll"<
exit(1);
}
result[1] = (f[0] * arr[2] - f[1] * arr[0])/det;
result[0] = (f[0]*arr[2] - arr[1]*arr[2]*result[1])/(arr[0]*arr[2]);
}
double norm22(double * f) {
return sqrt(f[0]* f[0] + f[1]* f[1]);
}
void computeG( double * g , double *x ,double a) {
g[0] = f1(x[0],x[1],a);
g[1] = f2(x[0],x[1],a);
}
void computeJ( double * j , double *x ,double a) {
//g[0] = f1(x[0],x[1],a);
//g[1] = f2(x[0],x[1],a);
j[0] = f1px1 (x[0],x[1],a);
j[1] = f1px2 (x[0],x[1],a);
j[2] = f2px1 (x[0],x[1],a);
j[3] = f2px2 (x[0],x[1],a);
}
bool isPrime( int n);
int main () {
double a = 173;
cout<<"a = "<
double g[2] = {0,0};
double j[4] = {0,0,0,0};
//double x[2] = {sqrt(a),sqrt(a)};
double x[2] = {sqrt(a),1.0};
double delta[2]= {0,0};
computeG( g, x,a);
computeJ( j, x,a);
//computeDeltaX (double * arr, double * f, double * result)
double ng[2] ={-1.0*g[0],-1.0*g[1]};
computeDeltaX (j, ng,delta) ;
cout<< norm22(delta) <
x[0] += delta[0];
x[1] += delta[1];
int i = 1;
while (! isVerySmall(norm22(delta)) && i<= __ITERATIONLIMIT) {
computeG( g, x,a);
computeJ( j, x,a);
ng[0] =-1.0*g[0];
ng[1] =-1.0*g[1];
//computeDeltaX (double * arr, double * f, double * result)
computeDeltaX (j, ng,delta) ;
x[0] += delta[0];
x[1] += delta[1];
cout<< norm22(delta) <
i++;
}
cout<<"result x : " <
cout<<"result g : " <
cout<<"iterations: "<
}
bool isPrime( int n) {
bool flag = true;
int l = (int)sqrt((double)n);
//cout<<"l:"<
int i = 0;
for(i = l; i>= 2;i-- ) {
if (n%i ==0 ) {
flag = false;
//cout<<"divisor:"<
break;
}
}
return flag;
} |
|
