#include <stdio.h>
#include <math.h>
#define DEFAULT_UPPER(10)
#define DEFAULT_LOWER(-10)
#define DEFAULT_E(0.00000001)
#define _MID(x,y)((x+y)/2)
#define _VALUE(x)(2*x*x*x-4*x*x+3*x-6)
double _e
int getRoot(double lower, double upper, double *result)
main()
{
double root
printf("Enter a deviation:")
scanf("%lf",&_e)
if(_e == 0.0)
_e = DEFAULT_E
if(getRoot(DEFAULT_LOWER, DEFAULT_UPPER, &root))
printf("Root:%2.8lf\n", root)
else
printf("Root:No Solution.\n")
}
int getRoot(double lower, double upper, double *result)
{
*result = _MID(lower,upper)
if(upper - lower <= _e)
return 1
if(_VALUE(lower)*_VALUE(*result) <= 0)
return getRoot(lower, *result, result)
else if(_VALUE(*result)*_VALUE(upper) <= 0)
return getRoot(*result, upper, result)
else
return 0
}
#include<iostream>#include<cstdlib>
#include<cmath>
using
namespace
std
#define
N
2
//用来设置方程组的行数
#define
eps
2.2204e-16
double*
MatrixMultiply(double*
J,double
Y[])
double*
Inv(double
*J)
double
norm(double
Q[])
double*
F(double
X[])
double*
JF(double
X[])
int
method(double*
Y,double
epsilon)
int
newdim(double
P[],double
delta,double
epsilon,int
max1,double
*err)
{
double
*Y=NULL,*J=NULL,Q[2]={0},*Z=NULL,*temp=NULL
double
relerr=0.0
int
k=0,i=0,iter=0
Y=F(P)
for(k=1k<max1k++)
{
J=JF(P)
temp=MatrixMultiply(Inv(J),Y)
for(i=0i<Ni++)
Q[i]=P[i]-temp[i]
Z=F(Q)
for(i=0i<Ni++)
temp[i]=Q[i]-P[i]
*err=norm(temp)
relerr=*err/(norm(Q)+eps)
for(i=0i<Ni++)
P[i]=Q[i]
for(i=0i<Ni++)
Y[i]=Z[i]
iter=k
if((*err<delta)||(relerr<delta)||method(Y,epsilon))
break
}
return
iter
}
int
method(double*
Y,double
epsilon)
{
if(fabs(Y[0])<epsilon&&fabs(Y[1])<epsilon)
return
1
else
return
0
}
//矩阵乘法，要求，J为方阵，Y为与J维数相同的列向量
double
*MatrixMultiply(double*
J,double
Y[])
{
double
*X=NULL
int
i=0,j=0
X=(double*)malloc(N*sizeof(double))
for(i=0i<Ni++)
X[i]=0
for(i=0i<Ni++)
for(j=0j<Nj++)
X[i]+=J[i*N+j]*Y[j]
return
X
}
//二阶矩阵的求逆（在M次多项式曲线拟合算法文件中给出了对任意可逆矩阵的求逆算法）
double
*Inv(double
*J)
{
double
X[4]={0},temp=0.0
int
i=0
temp=1/(J[0]*J[3]-J[1]*J[2])
X[0]=J[3]
X[1]=-J[1]
X[2]=-J[2]
X[3]=J[0]
for(i=0i<4i++)
J[i]=temp*X[i]
return
J
}
double
norm(double
Q[])
{
double
max=0.0
int
i=0
for(i=0i<Ni++)
{
if(Q[i]>max)
max=Q[i]
}
return
max
}
double*
F(double
X[])
{
double
x=X[0]
double
y=X[1]
double
*Z=NULL
Z=(double*)malloc(2*sizeof(double))
Z[0]=x*x-2*x-y+0.5
Z[1]=x*x+4*y*y-4
return
Z
}
double*
JF(double
X[])
{
double
x=X[0]
double
y=X[1]
double
*W=NULL
W=(double*)malloc(4*sizeof(double))
W[0]=2*x-2
W[1]=-1
W[2]=2*x
W[3]=8*y
return
W
}
main()
{
double
P[2]={0}
double
delta=0.0,epsilon=0.0,err=0.0
int
max1=0,iter=0,i=0
cout<<"牛顿法解非线性方程组:\nx^2-4-y+2=0\nx^2+4*y^2-2=0\n"
cout<<"\n输入的初始近似值x0,y0\n"
for(i=0i<2i++)
cin>>P[i]
cout<<"请依次输入P的误差限，F(P)的误差限，最大迭代次数\n"
cin>>delta
cin>>epsilon
cin>>err
iter=newdim(P,delta,epsilon,max1,&err)
cout<<"收敛到P的解为:\n"
for(i=0i<2i++)
cout<<"X("<<i+1<<")="<<P[i]<<endl
cout<<"\n迭代次数为:
"<<iter
cout<<"\n."<<err<<endl
getchar()
}
一种思路：将该方程组的系数组成一个矩阵，写成f(A)X=B的形式，然后可以像一元方程用Newton法来求解，只是注意导数是对矩阵的每一个元素求导，还要求逆矩阵
详细过程请找一本高等代数的教材就有了
