📄 gas_newton.cpp
字号:
#include "stdio.h"
#include "math.h"
#include "stdio.h"
#include "malloc.h"
#include <stdlib.h>
#include <string.h>
//void writeFile(double* data, int id);
double Surplus(double A[],int m,int n) /*行列式求值*/
{
int i,j,k,p,r;
double X,temp=1,temp1=1,s=0,s1=0;
if(n==2)
{
for(i=0;i<m;i++)
for(j=0;j<n;j++)
if((i+j)%2)temp1*=A[i*n+j];
else temp*=A[i*n+j];
X=temp-temp1;
}
else{
for(k=0;k<n;k++)
{
for(i=0,j=k;i<m,j<n;i++,j++)
temp*=A[i*n+j];
if(m-i)
{
for(p=m-i,r=m-1;p>0;p--,r--)
temp*=A[r*n+p-1];
}
s+=temp;
temp=1;
}
for(k=n-1;k>=0;k--)
{
for(i=0,j=k;i<m,j>=0;i++,j--)
temp1*=A[i*n+j];
if(m-i)
{
for(p=m-1,r=i;r<m;p--,r++)
temp1*=A[r*n+p];
}
s1+=temp1;
temp1=1;
}
X=s-s1;
}
return X;
}
int rinv(int n,double a[]) /* 矩阵求逆 */
{
double *is,*js;
int i,j,k,l,u,v;
double d,p;
is=(double*)malloc(n*sizeof(double));
js=(double*)malloc(n*sizeof(double));
for(k=0;k<=n-1;k++)
{
d=0.0;
for(i=k;i<=n-1;i++)
for(j=0;j<=n-1;j++)
{
l=i*n+j;
p=fabs(a[l]);
if(p>d){d=p;is[k]=i;js[k]=j;}
}
if(d+1.0==1.0)
{
free(is);
free(js);
printf("err**not inv\n");
return(0);
}
if(is[k]!=k)
for(j=0;j<=n-1;j++)/* 行交换*/
{
u=k*n+j;
v=is[k]*n+j;
p=a[u];
a[u]=a[v];
a[v]=p;
}
if(js[k]!=k)
for(i=0;i<=n-1;i++)/*列交换*/
{
u=i*n+k;
v=i*n+js[k];
p=a[u];
a[u]=a[v];
a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for(j=0;j<=n-1;j++)
if(j!=k)
{
u=k*n+j;
a[u]=a[u]*a[l];
}
for(i=0;i<=n-1;i++)
if(i!=k)
for(j=0;j<=n-1;j++)
if(j!=k)
{
u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for(i=0;i<=n-1;i++)
if(i!=k)
{
u=i*n+k;
a[u]=-a[u]*a[l];
};
}
for(k=n-1;k>=0;k--)
{
if(js[k]!=k)
for(j=0;j<=n-1;j++)
{
u=k*n+j;
v=js[k]*n+j;
p=a[u];
a[u]=a[v];
a[v]=p;
}
if(is[k]!=k)
for(i=0;i<=n-1;i++)
{
u=i*n+k;
v=i*n+is[k];
p=a[u];
a[u]=a[v];
a[v]=p;
}
}
free(is);
free(js);
return(1);
}
double *MatrixInver(double A[],int m,int n) /*矩阵转置*/
{
int i,j;
double *B=NULL;
B=(double *)malloc(m*n*sizeof(double));
for(i=0;i<n;i++)
for(j=0;j<m;j++)
B[i*m+j]=A[j*n+i];
return B;
}
double *fun(double M[],int n) /*通过该函数，计算归一化相速度及其偏导数并存放在一维数组Z[n+1]，返回数组首地址*/
{
double Q,A,B,C,Pxx,Pzz,Pxz;
double Pma1,Pma2,Pma3,Pma4,Pma5,Pma6,Pxxb1,Pxxb2,Pxxb3,Pxxc,Pxxd,Pzzb1,Pzzb2,Pzzb3,Pzzc,Pzzd,Pxze,Pxzf;
double Qa1,Qa2,Qa3,Qa4,Qa5,Qa6,Qb1,Qb2,Qb3,Qc,Qd,Qe,Qf;
double a1,a2,a3,a4,a5,a6,b1,b2,b3,c,d,e,f;
double a,b,kx,kz,D,Gs,pi=3.1415926; /*本行变量值已知*/
double T,T1,K1,K2,K3; /*中间变量*/
double *Z=NULL;
Z=(double*)malloc((n+1)*sizeof(double));
a=3000,b=1500,kx=kz=2*pi/100,D=10,Gs=5;
a1=M[0];
a2=M[1];
a3=M[2];
a4=M[3];
a5=M[4];
a6=M[5];
e=M[6];
f=M[7];
b1=M[8];
b2=M[9];
b3=M[10];
c=M[11];
d=M[12];
Pxx=-(4*c*pow(sin(kx*D/2),2)+d*pow(sin(kx*D),2))*(b1+2*b2*cos(kz*D)+2*b3*cos(2*kz*D));
Pzz=-(4*c*pow(sin(kz*D/2),2)+d*pow(sin(kz*D),2))*(b1+2*b2*cos(kx*D)+2*b3*cos(2*kx*D));
Pxz=-e*sin(kx*D)*sin(kz*D)-f/4*sin(2*kx*D)*sin(2*kz*D);
A=a1+2*a2*(cos(kx*D)+cos(kz*D))+4*a3*cos(kx*D)*cos(kz*D)+2*a4*(cos(2*kx*D)+cos(2*kz*D))+4*a5*(cos(2*kx*D)*cos(kz*D)+cos(kx*D)*cos(2*kz*D))+4*a6*cos(2*kx*D)*cos(2*kz*D);
B=-Pxx-Pzz;
C=pow((Pxx-Pzz),2)+4*pow(Pxz,2);
Pma1=1;
Pma2=2*(cos(kx*D)+cos(kz*D));
Pma3=4*cos(kx*D)*cos(kz*D);
Pma4=2*(cos(2*kx*D)+cos(2*kz*D));
Pma5=4*(cos(2*kx*D)*cos(kz*D)+cos(kx*D)*cos(2*kz*D));
Pma6=4*cos(2*kx*D)*cos(2*kz*D);
Pxze=-sin(kx*D)*sin(kz*D);
Pxzf=-sin(2*kx*D)*sin(2*kz*D)/4;
Pxxb1=-(4*c*pow(sin(kx*D/2),2)+d*pow(sin(kx*D),2));
Pxxb2=-(4*c*pow(sin(kx*D/2),2)+d*pow(sin(kx*D),2))*2*cos(kz*D);
Pxxb3=-(4*c*pow(sin(kx*D/2),2)+d*pow(sin(kx*D),2))*2*cos(2*kz*D);
Pxxc=-(b1+2*b2*cos(kz*D)+2*b3*cos(2*kz*D))*4*pow(sin(kx*D/2),2);
Pxxd=-(b1+2*b2*cos(kz*D)+2*b3*cos(2*kz*D))*pow(sin(kx*D),2);
Pzzb1=-(4*c*pow(sin(kz*D/2),2)+d*pow(sin(kz*D),2));
Pzzb2=-(4*c*pow(sin(kz*D/2),2)+d*pow(sin(kz*D),2))*2*cos(kx*D);
Pzzb3=-(4*c*pow(sin(kz*D/2),2)+d*pow(sin(kz*D),2))*2*cos(2*kx*D);
Pzzc=-(b1+2*b2*cos(kx*D)+2*b3*cos(2*kx*D))*4*pow(sin(kz*D/2),2);
Pzzd=-(b1+2*b2*cos(kx*D)+2*b3*cos(2*kx*D))*pow(sin(kz*D),2);
K1=1/(2*pi*b/a)*Gs; /*a,b分别为纵横波相速度，Gs为每个横波波长内的节点数*/
K2=1+b*b/(a*a);
K3=1-b*b/(a*a);
T=(1+b*b/(a*a))*B+(1-b*b/(a*a))*pow(C,1/2);
T1=1/(2*A)*T;
Q=K1*pow(T1,1/2); /*Q为纵波归一化相速度*/
Qa1=-K1/(4*A*A)*pow(T1,-1/2)*T*Pma1;
Qa2=-K1/(4*A*A)*pow(T1,-1/2)*T*Pma2;
Qa3=-K1/(4*A*A)*pow(T1,-1/2)*T*Pma3;
Qa4=-K1/(4*A*A)*pow(T1,-1/2)*T*Pma4;
Qa5=-K1/(4*A*A)*pow(T1,-1/2)*T*Pma5;
Qa6=-K1/(4*A*A)*pow(T1,-1/2)*T*Pma6;
Qe=K1/A*pow(T1,-1/2)*K3*Pxz*Pxze;
Qf=K1/A*pow(T1,-1/2)*K3*Pxz*Pxzf;
Qb1=K1/(4*A)*pow(T1,-1/2)*(-K2*(Pxxb1+Pzzb1)+K3*pow(C,-1/2)*(Pxx-Pzz)*(Pxxb1-Pzzb1));
Qb2=K1/(4*A)*pow(T1,-1/2)*(-K2*(Pxxb2+Pzzb2)+K3*pow(C,-1/2)*(Pxx-Pzz)*(Pxxb2-Pzzb2));
Qb3=K1/(4*A)*pow(T1,-1/2)*(-K2*(Pxxb3+Pzzb3)+K3*pow(C,-1/2)*(Pxx-Pzz)*(Pxxb3-Pzzb3));
Qc=K1/(4*A)*pow(T1,-1/2)*(-K2*(Pxxc+Pzzc)+K3*pow(C,-1/2)*(Pxx-Pzz)*(Pxxc-Pzzc));
Qd=K1/(4*A)*pow(T1,-1/2)*(-K2*(Pxxd+Pzzd)+K3*pow(C,-1/2)*(Pxx-Pzz)*(Pxxd-Pzzd));
Z[0]=Q;
Z[1]=Qa1;
Z[2]=Qa2;
Z[3]=Qa3;
Z[4]=Qa4;
Z[5]=Qa5;
Z[6]=Qa6;
Z[7]=Qe;
Z[8]=Qf;
Z[9]=Qb1;
Z[10]=Qb2;
Z[11]=Qb3;
Z[12]=Qc;
Z[13]=Qd;
return Z;
}
void main()
{
static double m1[13],m0[13]={0.51288,0.145,0.022,0.005,-0.00298,0.00011,1.20,
-0.026,0.61,0.27089,-0.02572,0.759,0.312};
double d,eps,p,dm[13]; /*m0初始化时依次用a1、a2、a3、a4、a5、a6、e、f、b1、b2、b3、c、d*/
double temp,sum,*p1;
int i,j,q,t,m=1,n=13;
double f0,jcb[13],e,jdm; /*jcb为m*n维矩阵*/
double arr1[169],arr2[13];
double arr7[13]={1,0,1,1,0,1,1,0,1,1,0,1,1};
static int num=0;
d=1.6;
eps=0.01;
e=eps+1.0;
while(e>=eps)
{
p1=fun(m0,n);
f0=*p1;
for(i=0;i<n;i++)
jcb[i]=*(p1+i+1); /*求得雅可比阵*/
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
sum=0.0;
for(t=0;t<m;t++)
{
temp=*(MatrixInver(jcb,m,n)+i*m+t)*jcb[t*n+j];
sum+=temp;
}
arr1[i*n+j]=sum; /*arr1为雅可比矩阵转置乘雅可比矩阵*/
}
}
q=rinv(n,arr1);
if(q!=0)
{
for(i=0;i<n;i++)
{
sum=0.0;
for(t=0;t<n;t++)
{
temp=*(arr1+i*n+t)**(MatrixInver(jcb,m,n)+t);
sum+=temp;
}
arr2[i]=sum; /*arr2为arr1的逆阵与雅可比阵转置之积*/
}
for(i=0;i<n;i++) /* 计算dm */
{
sum=0;
for(j=0;j<m;j++)
{
temp=arr2[i*m+j]*(d-f0);
sum+=temp;
}
dm[i]=sum;
}
sum=0.0; /* 计算jcb*dm */
for(i=0;i<n;i++)
{
temp=jcb[i]*dm[i];
sum+=temp;
}
jdm=sum;
for(i=0;i<n;i++) /* 计算m1 */
m1[i]=dm[i]+m0[i];
e=d-f0-jdm;
// if(e>p)p=e;
for(i=0;i<n;i++)
m0[i]=m1[i];
printf("f0=%12.7f\n",f0);
num++;
}
else break;
}
printf("a1=%12.7f\na2=%12.7f\na3=%12.7f\na4=%12.7f\na5=%12.7f\na6=%12.7f\ne=%12.7f\nf=%12.7f\nb1=%12.7f\nb2=%12.7f\nb3=%12.7f\nc=%12.7f\nd=%12.7f\n",m1[0],m1[1],m1[2],m1[3],m1[4],m1[5],m1[6],m1[7],m1[8],m1[9],m1[10],m1[11],m1[12]);
printf("\n");
printf("ee=%12.7f\n",e);
printf("num=%d\n",num);
printf("\n");
// writeFile(&m1[0], 13);
}
//void writeFile(double *data, int id)
//{
// FILE* file;
// file = fopen("readme2.txt","w");
//
//
// for(int i=0; i<id; i++)
// {
// char tmp[15] = {0};
//
// sprintf(tmp,"%12.7f",data[i]);
// fwrite(tmp, 15, 1, file);
// fwrite("\n", 1, 1, file);
// }
//
// fclose(file);
//}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -