matlab华军软件官方网站,用Forcal设计科学(计算)软件

#include

#include

#include

#include "math.h"

double k12,k21,k23,k32,k30;

double Array[3];

double aa[5]={1.0e-50,1.0e-50,1.0e-50,1.0e-50,1.0e-50};

double bb[5]={   1,      1,      1,      1,      1   };

double xx[6]={ 0.1,  0.1,  0.1,  0.1,  0.1};

double xxx[6][10];

double step, eps;

void gpbs1f(double t,double *y,int n,double *d)

{

d[0]=-k12*y[0]+k21*y[1];

d[1]=(-k21-k23)*y[1]+k12*y[0]+k32*y[2];

d[2]=(-k32-k30)*y[2]+k23*y[1];

return;

}

static void rkt1(double t,double h,int n,double *y,double *b,double *d,double *e)

{ int i,k;

double a[4],tt;

a[0]=h/2.0; a[1]=a[0]; a[2]=h; a[3]=h;

gpbs1f(t,y,n,d);

for (i=0; i<=n-1; i++) { b[i]=y[i]; e[i]=y[i];}

for (k=0; k<=2; k++)

{ for (i=0; i<=n-1; i++)

{ y[i]=e[i]+a[k]*d[i];

b[i]=b[i]+a[k+1]*d[i]/3.0;

}

tt=t+a[k];

gpbs1f(tt,y,n,d);

}

for (i=0; i<=n-1; i++)

y[i]=b[i]+h*d[i]/6.0;

return;

}

void gpbs1(double t,double h,int n,double *y,double eps)

{ int i,j,k,m,nn,it;

double x,hh,dd,q,p,g[10],*b,*d,*u,*v,*w,*e;

b=(double *)malloc(10*n*sizeof(double));

d=(double *)malloc(n*sizeof(double));

u=(double *)malloc(n*sizeof(double));

v=(double *)malloc(n*sizeof(double));

w=(double *)malloc(n*sizeof(double));

e=(double *)malloc(n*sizeof(double));

for (j=0; j<=n-1; j++) v[j]=y[j];

x=t; nn=1; hh=h; g[0]=hh;

rkt1(x,hh,n,y,w,d,e);

for (j=0; j<=n-1; j++)

{ b[j]=y[j]; u[j]=y[j];}

k=1; it=1;

while (it==1)

{ nn=nn+nn; hh=hh/2.0; it=0;

g[k]=hh;

for (j=0; j<=n-1; j++) y[j]=v[j];

t=x;

for (j=0; j<=nn-1; j++)

{ rkt1(t,hh,n,y,w,d,e);

t=t+hh;

}

for (j=0; j<=n-1; j++)

{ dd=y[j]; m=0;

for (i=0; i<=k-1; i++)

if (m==0)

{ q=dd-b[i*n+j];

if (fabs(q)+1.0==1.0) m=1;

else dd=(g[k]-g[i])/q;

}

b[k*n+j]=dd;

if (m!=0) b[k*n+j]=1.0e+35;

}

for (j=0; j<=n-1; j++)

{ dd=0.0;

for (i=k-1; i>=0; i--)

dd=-g[i]/(b[(i+1)*n+j]+dd);

y[j]=dd+b[j];

}

p=0.0;

for (j=0; j<=n-1; j++)

{ q=fabs(y[j]-u[j]);

if (q>p) p=q;

}

if ((p>=eps)&&(k<7))

{ for (j=0; j<=n-1; j++) u[j]=y[j];

k=k+1; it=1;

}

}

free(b); free(d); free(u); free(v); free(w); free(e);

return;

}

double rn(double *rr)

{ int m;

double y,s,u,v;

s=65536.0; u=2053.0; v=13849.0;

*rr=u*(*rr)+v; m=(int)(*rr/s); *rr=*rr-m*s;

y=*rr/s;

return(y);

}

void t_i(double *a,double step,double eps,double &t1,double t2,double &x1,double &x2,double &x3) //用于约束条件定义

{   double h;

h=(t2-t1)/step;

do{   gpbs1(t1,h,3,a,eps);  //连分式法对微分方程组积分一步函数gpbs1

t1=t1+h;

}while(abs(t1-t2)>=h/2);

x1=a[0];x2=a[1];x3=a[2];

};

double t_i_2(double *a,double step,double eps,double &t1,double t2,double x_1,double x_2,double x_3) //用于计算目标函数

{   double h;

h=(t2-t1)/step;

do{   gpbs1(t1,h,3,a,eps);  //连分式法对微分方程组积分一步函数gpbs1

t1=t1+h;

}while(abs(t1-t2)>=h/2);

return (a[0]-x_1)*(a[0]-x_1)+(a[1]-x_2)*(a[1]-x_2)+(a[2]-x_3)*(a[2]-x_3);

};

double jcplxf(double *x,int n)

{ double t1,s=0;

k12=x[0];k21=x[1];k23=x[2];k32=x[3];k30=x[4];

t1=0; Array[0]=100; Array[1]=0; Array[2]=0;

s=s+ t_i_2(Array,step,eps,t1, 1 , 90,  8,  2);

s=s+ t_i_2(Array,step,eps,t1, 3 , 73, 19,  7);

return s+t_i_2(Array,step,eps,t1, 5 , 60, 14, 23);

}

void jcplxs(int n,int m,double *x,double *c,double *d,double *w)

{ double x1,x2,x3,t1;

c[0]=0.0; c[1]=0.0; c[2]=0.0;

d[0]=100; d[1]=100; d[2]=100;

k12=x[0];k21=x[1];k23=x[2];k32=x[3];k30=x[4];

t1=0; Array[0]=100; Array[1]=0; Array[2]=0;

t_i(Array,step,eps,t1, 1 ,x1, x2, x3); w[0]=x1+x2+x3;

t_i(Array,step,eps,t1, 3 ,x1, x2, x3); w[1]=x1+x2+x3;

t_i(Array,step,eps,t1, 5 ,x1, x2, x3); w[2]=x1+x2+x3;

return;

}

int jcplx(int n,int m,double *a,double *b,double alpha,double eps,double *x,double *xx,int k,double dx)

{

int r,g,i,j,it,kt,jt,kk;

double fj,fr,fg,z,rr,*c,*d,*w,*xt,*xf;

c=(double *)malloc(m*sizeof(double));

d=(double *)malloc(m*sizeof(double));

w=(double *)malloc(m*sizeof(double));

xt=(double *)malloc(n*sizeof(double));

xf=(double *)malloc(n*sizeof(double));

rr=0.0;

for (i=0; i<=n-1; i++)

xx[i*n*2]=x[i];

xx[n*n*2]=jcplxf(x,n);

for (j=1; j<=2*n-1; j++)

{ for (i=0; i<=n-1; i++)

{ xx[i*n*2+j]=a[i]+(b[i]-a[i])*rn(&rr);

x[i]=xx[i*n*2+j];

}

it=1;

while (it==1)

{ it=0; r=0; g=0;

while ((r

{ if ((a[r]<=x[r])&&(b[r]>=x[r])) r=r+1;

else g=1;

}

if (g==0)

{ jcplxs(n,m,x,c,d,w);

r=0;

while ((r

{ if ((c[r]<=w[r])&&(d[r]>=w[r])) r=r+1;

else g=1;

}

}

if (g!=0)

{ for (r=0; r<=n-1; r++)

{ z=0.0;

for (g=0; g<=j-1; g++)

z=z+xx[r*n*2+g]/(1.0*j);

xx[r*n*2+j]=(xx[r*n*2+j]+z)/2.0;

x[r]=xx[r*n*2+j];

}

it=1;

}

else xx[n*n*2+j]=jcplxf(x,n);

}

}

kk=1; it=1;

while (it==1)

{ it=0;

fr=xx[n*n*2]; r=0;

for (i=1; i<=2*n-1; i++)

if (xx[n*n*2+i]>fr)

{ r=i; fr=xx[n*n*2+i];}

g=0; j=0; fg=xx[n*n*2];

if (r==0)

{ g=1; j=1; fg=xx[n*n*2+1];}

for (i=j+1; i<=2*n-1; i++)

if (i!=r)

if (xx[n*n*2+i]>fg)

{ g=i; fg=xx[n*n*2+i];}

for (i=0; i<=n-1; i++)

{ xf[i]=0.0;

for (j=0; j<=2*n-1; j++)

if (j!=r)

xf[i]=xf[i]+xx[i*n*2+j]/(2.0*n-1.0);

xt[i]=(1.0+alpha)*xf[i]-alpha*xx[i*n*2+r];

}

jt=1;

while (jt==1)

{ jt=0;

z=jcplxf(xt,n);

while (z>fg)

{ for (i=0; i<=n-1; i++)

xt[i]=(xt[i]+xf[i])/2.0;

z=jcplxf(xt,n);

}

j=0;

for (i=0; i<=n-1; i++)

{ if (a[i]>xt[i])

{ xt[i]=xt[i]+dx; j=1;}

if (b[i]

{ xt[i]=xt[i]-dx; j=1;}

}

if (j!=0) jt=1;

else

{ jcplxs(n,m,xt,c,d,w);

j=0; kt=1;

while ((kt==1)&&(j

{ if ((c[j]<=w[j])&&(d[j]>=w[j])) j=j+1;

else kt=0;

}

if (j

{ for (i=0; i<=n-1; i++)

xt[i]=(xt[i]+xf[i])/2.0;

jt=1;

}

}

}

for (i=0; i<=n-1; i++)

xx[i*n*2+r]=xt[i];

xx[n*n*2+r]=z;

fr=0.0; fg=0.0;

for (j=0; j<=2*n-1; j++)

{ fj=xx[n*n*2+j];

fr=fr+fj/(2.0*n);

fg=fg+fj*fj;

}

fr=(fg-2.0*n*fr*fr)/(2.0*n-1.0);

if (fr>=eps)

{ kk=kk+1;

if (kk

}

}

for (i=0; i<=n-1; i++)

{ x[i]=0.0;

for (j=0; j<=2*n-1; j++)

x[i]=x[i]+xx[i*n*2+j]/(2.0*n);

}

z=jcplxf(x,n); x[n]=z;

free(c); free(d); free(w);

free(xt); free(xf);

return(kk);

}

int main(int argc, char *argv[])

{

clock_t tm;

double _eps,alpha,dx;

int k,i;

tm=clock();

step=20;eps=1e-6; //step越大，eps越小越精确，用于积分一步的变步长欧拉方法

_eps=1e-10; alpha=1.01;k=800;dx=1e-4;   //_eps控制精度要求；alpha为反射系数；k为允许的最大迭代次数；dx为微小增量。需选择合适的alpha,_eps,k,dx值求解。

i=jcplx(5,3,aa,bb,alpha,_eps,xx,&xxx[0][0],k,dx);

printf("\r\n实际迭代次数=%d\r\n",i);

printf("%12.6e, %12.6e, %12.6e, %12.6e, %12.6e, %12.6e\n",xx[0],xx[1],xx[2],xx[3],xx[4],xx[5]);

printf("程序运行 %e 秒。\n", double(clock()-tm)/1000.0);

}