#include "stdafx.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "NonLinearEquation.h"
#include "MatrixAlgo.h"
#include "LinearEquation.h"
#include "MatrixEigenvalue.h"
//求非线性方程一个实根的连分式解法
int dpqrt(double *x, double eps)
{
int i,j,m,it,l;
double a[10],y[10],z,h,x0,q;
l=10; q=1.0e+35; x0=*x; h=0.0;
while (l!=0)
{
l=l-1;
j=0;
it=l;
while (j<=7)
{
if (j<=2)
z=x0+0.1*j;
else
z=h;
y[j]=dpqrtf(z);
h=z;
if (j==0)
a[0]=z;
else
{
m=0;
i=0;
while ((m==0)&&(i<=j-1))
{
if (fabs(h-a[i])+1.0==1.0)
m=1;
else
h=(y[j]-y[i])/(h-a[i]);
i=i+1;
}
a[j]=h;
if (m!=0)
a[j]=q;
h=0.0;
for (i=j-1; i>=0; i--)
{
if (fabs(a[i+1]+h)+1.0==1.0)
h=q;
else
h=-y[i]/(a[i+1]+h);
}
h=h+a[0];
}
if (fabs(y[j])>=eps)
j=j+1;
else
{
j=10;
l=0;
}
}
x0=h;
}
*x=h;
return(10-it);
}
//求非线性方程一个实根的牛顿法
int dnewt(double *x, double eps, int js)
{
int k,l;
double y[2],d,p,x0,x1;
l=js;
k=1;
x0=*x;
dnewtf(x0,y);
d=eps+1.0;
while ((d>=eps)&&(l!=0))
{
if (fabs(y[1])+1.0==1.0)
{
printf("err\n");
return(-1);
}
x1=x0-y[0]/y[1];
dnewtf(x1,y);
d=fabs(x1-x0);
p=fabs(y[0]);
if (p>d)
d=p;
x0=x1;
l=l-1;
}
*x=x1;
k=js-l;
return(k);
}
//求复系数代数方程全部根的牛顿-下山法
int dcsrt(double ar[], double ai[], int n, double xr[], double xi[])
{
int m,i,jt,k,is,it;
double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
double g,u,v,pq,g1,u1,v1;
m=n;
p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
while ((m>0)&&(p+1.0==1.0))
{
m=m-1;
p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
}
if (m<=0)
{
printf("fail\n");
return(-1);
}
for (i=0; i<=m; i++)
{
ar[i]=ar[i]/p;
ai[i]=ai[i]/p;
}
for (i=0; i<=m/2; i++)
{
w=ar[i];
ar[i]=ar[m-i];
ar[m-i]=w;
w=ai[i];
ai[i]=ai[m-i];
ai[m-i]=w;
}
k=m; is=0; w=1.0;
jt=1;
while (jt==1)
{
pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
while (pq<1.0e-12)
{
xr[k-1]=0.0;
xi[k-1]=0.0;
k=k-1;
if (k==1)
{
p=ar[0]*ar[0]+ai[0]*ai[0];
xr[0]=-w*(ar[0]*ar[1]+ai[0]*ai[1])/p;
xi[0]=w*(ar[1]*ai[0]-ar[0]*ai[1])/p;
return(1);
}
pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
}
q=log(pq);
q=q/(1.0*k);
q=exp(q);
p=q;
w=w*p;
for (i=1; i<=k; i++)
{
ar[i]=ar[i]/q;
ai[i]=ai[i]/q;
q=q*p;
}
x=0.0001;
x1=x;
y=0.2;
y1=y;
dx=1.0;
g=1.0e+37;
l40:
u=ar[0];
v=ai[0];
for (i=1; i<=k; i++)
{
p=u*x1;
q=v*y1;
pq=(u+v)*(x1+y1);
u=p-q+ar[i];
v=pq-p-q+ai[i];
}
g1=u*u+v*v;
if (g1>=g)
{
if (is!=0)
{
it=1;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
else
{
g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
if (t>=1.0e-03)
goto l40;
if (g>1.0e-18)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
}
g90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
}
else
{
g=g1;
x=x1;
y=y1;
is=0;
if (g<=1.0e-22)
g90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
else
{
u1=k*ar[0];
v1=ai[0];
for (i=2; i<=k; i++)
{
p=u1*x;
q=v1*y;
pq=(u1+v1)*(x+y);
u1=p-q+(k-i+1)*ar[i-1];
v1=pq-p-q+(k-i+1)*ai[i-1];
}
p=u1*u1+v1*v1;
if (p<=1.0e-20)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,
&c,&k,&is,&it);
if (it==0)
goto l40;
g90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
}
else
{
dx=(u*u1+v*v1)/p;
dy=(u1*v-v1*u)/p;
t=1.0+4.0/k;
g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
if (t>=1.0e-03)
goto l40;
if (g>1.0e-18)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
g90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
}
}
}
if (k==1)
jt=0;
else
jt=1;
}
return(1);
}
void g60(double *t, double *x, double *y, double *x1, double *y1, double *dx, double *dy, double *p, double *q, int *k, int *it)
{
*it=1;
while (*it==1)
{
*t=*t/1.67;
*it=0;
*x1=*x-*t*(*dx);
*y1=*y-*t*(*dy);
if (*k>=30)
{
*p=sqrt(*x1*(*x1)+*y1*(*y1));
*q=exp(75.0/(*k));
if (*p>=*q)
*it=1;
}
}
return;
}
void g90(double xr[], double xi[], double ar[], double ai[], double *x, double *y, double *p, double *w, int *k)
{
int i;
for (i=1; i<=*k; i++)
{
ar[i]=ar[i]+ar[i-1]*(*x)-ai[i-1]*(*y);
ai[i]=ai[i]+ar[i-1]*(*y)+ai[i-1]*(*x);
}
xr[*k-1]=*x*(*w);
xi[*k-1]=*y*(*w);
*k=*k-1;
if (*k==1)
{
*p=ar[0]*ar[0]+ai[0]*ai[0];
xr[0]=-*w*(ar[0]*ar[1]+ai[0]*ai[1])/(*p);
xi[0]=*w*(ar[1]*ai[0]-ar[0]*ai[1])/(*p);
}
return;
}
void g65(double *x, double *y, double *x1, double *y1, double *dx, double *dy, double *dd, double *dc, double *c, int *k, int *is, int *it)
{
if (*it==0)
{
*is=1;
*dd=sqrt(*dx*(*dx)+*dy*(*dy));
if (*dd>1.0)
*dd=1.0;
*dc=6.28/(4.5*(*k));
*c=0.0;
}
while(1==1)
{
*c=*c+*dc;
*dx=*dd*cos(*c);
*dy=*dd*sin(*c);
*x1=*x+*dx;
*y1=*y+*dy;
if (*c<=6.29)
{
*it=0;
return;
}
*dd=*dd/1.67;
if (*dd<=1.0e-07)
{
*it=1;
return;
}
*c=0.0;
}
}
//求实系数代数方程全部报的牛顿-下山法
int dsrrt(double a[], int n, double xr[], double xi[])
{
int m,i,jt,k,is,it;
double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
double g,u,v,pq,g1,u1,v1;
m=n;
while ((m>0)&&(fabs(a[m])+1.0==1.0))
m=m-1;
if (m<=0)
{
printf("fail\n");
return(-1);
}
for (i=0; i<=m; i++)
a[i]=a[i]/a[m];
for (i=0; i<=m/2; i++)
{
w=a[i];
a[i]=a[m-i];
a[m-i]=w;
}
k=m;
is=0;
w=1.0;
jt=1;
while (jt==1)
{
pq=fabs(a[k]);
while (pq<1.0e-12)
{
xr[k-1]=0.0;
xi[k-1]=0.0;
k=k-1;
if (k==1)
{
xr[0]=-a[1]*w/a[0];
xi[0]=0.0;
return(1);
}
pq=fabs(a[k]);
}
q=log(pq);
q=q/(1.0*k);
q=exp(q);
p=q; w=w*p;
for (i=1; i<=k; i++)
{
a[i]=a[i]/q;
q=q*p;
}
x=0.0001;
x1=x;
y=0.2;
y1=y;
dx=1.0;
g=1.0e+37;
l40:
u=a[0];
v=0.0;
for (i=1; i<=k; i++)
{
p=u*x1; q=v*y1;
pq=(u+v)*(x1+y1);
u=p-q+a[i];
v=pq-p-q;
}
g1=u*u+v*v;
if (g1>=g)
{
if (is!=0)
{
it=1;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
else
{
g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
if (t>=1.0e-03)
goto l40;
if (g>1.0e-18)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
}
g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
}
else
{
g=g1;
x=x1;
y=y1;
is=0;
if (g<=1.0e-22)
g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
else
{
u1=k*a[0];
v1=0.0;
for (i=2; i<=k; i++)
{
p=u1*x;
q=v1*y;
pq=(u1+v1)*(x+y);
u1=p-q+(k-i+1)*a[i-1];
v1=pq-p-q;
}
p=u1*u1+v1*v1;
if (p<=1.0e-20)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
}
else
{
dx=(u*u1+v*v1)/p;
dy=(u1*v-v1*u)/p;
t=1.0+4.0/k;
g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
if (t>=1.0e-03)
goto l40;
if (g>1.0e-18)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
}
}
}
if (k==1)
jt=0;
else
jt=1;
}
return(1);
}
//求非线性方程组一组实根的蒙特卡络法
void dnmtc(double x[], int n, double b, int m, double eps)
{
int k,i;
double a,r,*y,z,z1;
y=(double *)malloc(n*sizeof(double));
a=b;
k=1;
r=1.0;
z=dnmtcf(x,n);
while (a>=eps)
{
for (i=0; i<=n-1; i++)
y[i]=-a+2.0*a*mrnd1(&r)+x[i];
z1=dnmtcf(y,n);
k=k+1;
if (z1>=z)
{
if (k>m)
{
k=1;
a=a/2.0;
}
}
else
{
k=1;
for (i=0; i<=n-1; i++)
x[i]=y[i];
z=z1;
if (z<eps)
return;
}
}
return;
}
//求非线性方程组一组实根的拟牛顿法
int dnetn(int n, double eps, double t, double h, double x[], int k)
{
int i,j,l;
double am,z,beta,d,*y,*a,*b;
y=(double *)malloc(n*sizeof(double));
a=(double *)malloc(n*n*sizeof(double));
b=(double *)malloc(n*sizeof(double));
l=k; am=1.0+eps;
while (am>=eps)
{
dnetnf(x,b,n);
am=0.0;
for (i=0; i<=n-1; i++)
{
z=fabs(b[i]);
if (z>am) am=z;
}
if (am>=eps)
{
l=l-1;
if (l==0)
{
free(y);
free(b);
free(a);
printf("fail\n");
return(0);
}
for (j=0; j<=n-1; j++)
{
z=x[j];
x[j]=x[j]+h;
dnetnf(x,y,n);
for (i=0; i<=n-1; i++)
a[i*n+j]=y[i];
x[j]=z;
}
if (agaus(a,b,n)==0)
{
free(y);
free(a);
free(b);
return(-1);
}
beta=1.0;
for (i=0; i<=n-1; i++)
beta=beta-b[i];
if (fabs(beta)+1.0==1.0)
{
free(y);
free(a);
free(b);
printf("fail\n");
return(-2);
}
d=h/beta;
for (i=0; i<=n-1; i++)
x[i]=x[i]-d*b[i];
h=t*h;
}
}
free(y);
free(a);
free(b);
return(k-l);
}
//求非线性方程组一组实根的梯度法
int dsnse(int n, double eps, double x[], int js)
{
int l,j;
double f,d,s,*y;
y=(double *)malloc(n*sizeof(double));
l=js;
f=dsnsef(x,y,n);
while (f>=eps)
{
l=l-1;
if (l==0)
{
free(y);
return(js);
}
d=0.0;
for (j=0; j<=n-1; j++)
d=d+y[j]*y[j];
if (d+1.0==1.0)
{
free(y);
return(-1);
}
s=f/d;
for (j=0; j<=n-1; j++)
x[j]=x[j]-s*y[j];
f=dsnsef(x,y,n);
}
free(y);
return(js-l);
}
//求非线性方程一个实根的蒙特卡洛法
void dmtcl(double *x, double b, int m, double eps)
{
int k;
double xx,a,r,y,x1,y1;
a=b;
k=1;
r=1.0;
xx=*x;
y=dmtclf(xx);
while (a>=eps)
{
x1=mrnd1(&r);
x1=-a+2.0*a*x1;
x1=xx+x1;
y1=dmtclf(x1);
k=k+1;
if (fabs(y1)>=fabs(y))
{
if (k>m)
{
k=1;
a=a/2.0;
}
}
else
{
k=1;
xx=x1;
y=y1;
if (fabs(y)<eps)
{
*x=xx;
return;
}
}
}
*x=xx;
return;
}
//求实系数代数方程全部根的QR方法
int dqrrt(double a[], int n, double xr[], double xi[], double eps, int jt)
{
int i,j;
double *q;
q=(double *)malloc(n*n*sizeof(double));
for (j=0; j<=n-1; j++)
q[j]=-a[n-j-1]/a[n];
for (j=n; j<=n*n-1; j++)
q[j]=0.0;
for (i=0; i<=n-2; i++)
q[(i+1)*n+i]=1.0;
i=chhqr(q,n,xr,xi,eps,jt);
free(q);
return(i);
}
//求非线性方程组最小二乘解的广义逆法
int dngin(int m, int n, double eps1, double eps2, double x[], int ka)
{
int i,j,k,l,kk,jt;
double y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
double *p,*d,*pp,*dx,*u,*v,*w;
p=(double *)malloc(m*n*sizeof(double));
d=(double *)malloc(m*sizeof(double));
pp=(double *)malloc(n*m*sizeof(double));
dx=(double *)malloc(n*sizeof(double));
u=(double *)malloc(m*m*sizeof(double));
v=(double *)malloc(n*n*sizeof(double));
w=(double *)malloc(ka*sizeof(double));
l=60; alpha=1.0;
while (l>0)
{
dnginf(m,n,x,d);
dngins(m,n,x,p);
jt=agmiv(p,m,n,d,dx,pp,eps2,u,v,ka);
if (jt<0)
{
free(p);
free(d);
free(pp);
free(w);
free(dx);
free(u);
free(v);
return(jt);
}
j=0;
jt=1;
h2=0.0;
while (jt==1)
{
jt=0;
if (j<=2)
z=alpha+0.01*j;
else
z=h2;
for (i=0; i<=n-1; i++)
w[i]=x[i]-z*dx[i];
dnginf(m,n,w,d);
y1=0.0;
for (i=0; i<=m-1; i++)
y1=y1+d[i]*d[i];
for (i=0; i<=n-1; i++)
w[i]=x[i]-(z+0.00001)*dx[i];
dnginf(m,n,w,d);
y2=0.0;
for (i=0; i<=m-1; i++)
y2=y2+d[i]*d[i];
y0=(y2-y1)/0.00001;
if (fabs(y0)>1.0e-10)
{
h1=y0;
h2=z;
if (j==0)
{
y[0]=h1;
b[0]=h2;
}
else
{
y[j]=h1;
kk=0;
k=0;
while ((kk==0)&&(k<=j-1))
{
y3=h2-b[k];
if (fabs(y3)+1.0==1.0)
kk=1;
else
h2=(h1-y[k])/y3;
k=k+1;
}
b[j]=h2;
if (kk!=0)
b[j]=1.0e+35;
h2=0.0;
for (k=j-1; k>=0; k--)
h2=-y[k]/(b[k+1]+h2);
h2=h2+b[0];
}
j=j+1;
if (j<=7)
jt=1;
else
z=h2;
}
}
alpha=z;
y1=0.0;
y2=0.0;
for (i=0; i<=n-1; i++)
{
dx[i]=-alpha*dx[i];
x[i]=x[i]+dx[i];
y1=y1+fabs(dx[i]);
y2=y2+fabs(x[i]);
}
if (y1<eps1*y2)
{
free(p);
free(pp);
free(d);
free(w);
free(dx);
free(u);
free(v);
return(1);
}
l=l-1;
}
free(p);
free(pp);
free(d);
free(dx);
free(u);
free(v);
free(w);
return(0);
}
//求非线性方程实根的对分法
int ddhrt(double a, double b, double h, double eps, double x[], int m)
{
int n,js;
double z,y,z1,y1,z0,y0;
n=0;
z=a;
y=ddhrtf(z);
while ((z<=b+h/2.0)&&(n!=m))
{
if (fabs(y)<eps)
{
n=n+1;
x[n-1]=z;
z=z+h/2.0;
y=ddhrtf(z);
}
else
{
z1=z+h;
y1=ddhrtf(z1);
if (fabs(y1)<eps)
{
n=n+1;
x[n-1]=z1;
z=z1+h/2.0;
y=ddhrtf(z);
}
else if (y*y1>0.0)
{
y=y1;
z=z1;
}
else
{
js=0;
while (js==0)
{
if (fabs(z1-z)<eps)
{
n=n+1;
x[n-1]=(z1+z)/2.0;
z=z1+h/2.0;
y=ddhrtf(z);
js=1;
}
else
{
z0=(z1+z)/2.0;
y0=ddhrtf(z0);
if (fabs(y0)<eps)
{
x[n]=z0;
n=n+1;
js=1;
z=z0+h/2.0;
y=ddhrtf(z);
}
else if ((y*y0)<0.0)
{
z1=z0;
y1=y0;
}
else
{
z=z0;
y=y0;
}
}
}
}
}
}
return(n);
}
//求实函数或复函数方程一个复根的蒙特卡洛法
void dcmtc(double *x, double *y, double b, int m, double eps)
{
int k;
double xx,yy,a,r,z,x1,y1,z1;
a=b;
k=1;
r=1.0;
xx=*x;
yy=*y;
z=dcmtcf(xx,yy);
while (a>=eps)
{
x1=-a+2.0*a*mrnd1(&r);
x1=xx+x1;
y1=-a+2.0*a*mrnd1(&r);
y1=yy+y1;
z1=dcmtcf(x1,y1);
k=k+1;
if (z1>=z)
{
if (k>m)
{
k=1;
a=a/2.0;
}
}
else
{
k=1;
xx=x1;
yy=y1;
z=z1;
if (z<eps)
{
*x=xx;
*y=yy;
return;
}
}
}
*x=xx;
*y=yy;
return;
}
//求非线性方程一个实根的埃待金达代法
int datkn(double *x, double eps, int js)
{
int flag,l;
double u,v,x0;
l=0;
x0=*x;
flag=0;
while ((flag==0)&&(l!=js))
{
l=l+1;
u=datknf(x0);
v=datknf(u);
if (fabs(u-v)<eps)
{
x0=v;
flag=1;
}
else
x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
}
*x=x0;
l=js-l;
return(l);
}
double ddhrtf(double x)
{
double z;
z = (((((x-5.0)*x + 3.0)*x + 1.0) * x - 7.0) * x + 7.0) * x - 20.0;;
return z;
}
void dnewtf(double x, double y[2])
{
y[0] = x * x * (x - 1.0) - 1.0;
y[1] = 3.0 * x * x - 2.0 * x;
}
double datknf(double x)
{
double z;
z = 6.0 - x * x;
return z;
}
double dpqrtf(double x)
{
double z;
z = x * x * (x - 1.0) - 1.0;
return z;
}
/*
f[0] = x[0] - 5 * x[1]^2 + 7 * x[2]^2 + 12 = 0
f[1] = 3 * x[0] * x[1] + x[0] * x[2] - 11 * x[0] = 0
f[2] = 2 * x[1] * x[2] + 40 * x[0] = 0
*/
double dsnsef(double x[], double y[], int n)
{
double z, f1, f2, f3, df1, df2, df3;
n = n;
f1 = x[0] - 5.0 * x[1] * x[1] + 7.0 * x[2] + 12.0;
f2 = 3.0 * x[0] * x[1] + x[0] * x[2] - 11.0 * x[0];
f3 = 2.0 * x[1] * x[2] + 40.0 * x[0];
z = f1 * f1 + f2 * f2 + f3 * f3;
df1 = 1.0;
df2 = 3.0 * x[1] + x[2] - 11.0;
df3 = 40.0;
y[0] = 2.0 * (f1 * df1 + f2 * df2 + f3 * df3);
df1 = 10.0 * x[1];
df2 = 3.0 * x[0];
df3 = 2.0 * x[2];
y[1] = 2.0 * (f1 * df1 + f2 * df2 + f3 * df3);
df1 = 14.0 * x[2];
df2 = x[0];
df3 = 2.0 * x[1];
y[2] = 2.0 * (f1 * df1 + f2 * df2 + f3 * df3);
return z;
}
/*
f[0] = x[0] ^ 2 + x[1] ^ 2 + x[2] ^ 2 - 1 = 0
f[1] = 2 * x[0] ^ 2 + x[1] ^ 2 - 4 * x[2] = 0
f[3] = 3 * x[0] ^ 2 - 4 * x[1] + x[2] ^ 2 = 0
*/
void dnetnf(double x[], double y[], int n)
{
y[0] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] - 1;
y[1] = 2 * x[0] * x[0] + x[1] * x[1] - 4 * x[2];
y[3] = 3 * x[0] * x[0] - 4 * x[1] + x[2] * x[2];
}
/*
f[0] = x[0] ^ 2 + 10 * x[0] * x[1] + 4 * x[1] ^ 2 + 0.7401006 = 0
f[1] = x[0] ^ 2 - 3 * x[0] * x[1] + 2 * x[1] - 1.0201228 ＝ 0
其中m = 2, n = 2。取初值为（0.5，-1.0），eps1 = 0.000001，eps2 = 0.000001。
*/
void dnginf(int m, int n, double x[], double d[])
{
m = m;
n = n;
d[0] = x[0] * x[0] + 10 * x[0] * x[1] + 4 * x[1] * x[1] + 0.7401006;
d[1] = x[0] * x[0] - 3 * x[0] * x[1] + 2 * x[1] - 1.0201228;
}
void dngins(int m, int n, double x[], double p[])
{
m = m;
n = n;
p[0] = 2.0 * x[0] + 10.0 * x[1];
p[1] = 10.0 * x[0] + 8.0 * x[1];
p[2] = 2.0 * x[0] - 3.0 * x[1];
p[3] = -3.0 * x[0] + 4.0 * x[1];
}
/*
f(x) = e ^ (-x ^ 3) - sin(x) / cos(x) + 800 = 0
*/
double dmtclf(double x)
{
double z;
z = exp(-x * x * x) - sin(x)/cos(x) + 800.0;
return z;
}
/*
f(z) ＝z ^ 2 - 6z + 13 ＝ 0
取初值z = 0.5 + 0.5j, b = 1.0, m = 10, ε= 0.000001。其中
f(z) = (x ^ 2 - y ^ 2 - 6x + 13) + j(2xy - 6y)
*/
double dcmtcf(double x, double y)
{
double u, v, z;
u = x * x - y * y - 6.0 * x + 13.0;
v = 2.0 * x * y - 6.0 * y;
z = sqrt(u * u + v * v);
return z;
}
/*
f[0] = 3x[1] + x[2] + 2x[3] ^ 2 - 3 = 0;
f[1] = -3x[1] + 5x[2] ^ 2 + 2x[1]x[3] - 1 = 0;
f[2] = 25x[1]x[2] + 20x[3] + 12 = 0;
取初值X = (0, 0, 0)^T, b = 2.0, m = 10, ε = 0.00001。
*/
double dnmtcf(double x[], int n)
{
double f, f1, f2, f3;
n = n;
f1 = 3.0 * x[0] + x[1] + 2.0 * x[2] * x[2] - 3.0;
f2 = -3.0 * x[0] + 5.0 * x[1] * x[1] + 2.0 * x[0] * x[2] - 1.0;
f3 = 25.0 * x[0] * x[1] + 20.0 * x[2] + 12.0;
f = sqrt(f1 * f1 + f2 * f2 + f3 * f3);
return f;
}
----根据《C语言常用算法程序集》整理
扫一扫
491
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
