📄 type.txt
字号:
for(i=k+1;i<=n;i++) l[i-1][k-1]=s[i-1]/u[k-1][k-1];
}
for(k=1;k<=n-1;k++) {
temp=b[k-1];
b[k-1]=b[m[k-1]-1];
b[m[k-1]-1]=temp;
}
y[0]=b[0];
for(i=2;i<=n;i++){
l_u=0;
for(t=1;t<=i-1;t++) l_u=l_u+l[i-1][t-1]*y[t-1];
y[i-1]=b[i-1]-l_u;
}
x[n-1]=y[n-1]/u[n-1][n-1];
for(i=n-1;i>=1;i--){
u_x=0;
for(t=i+1;t<=n;t++)u_x=u_x+u[i-1][t-1]*x[t-1];
x[i-1]=(y[i-1]-u_x)/u[i-1][i-1];
}
do{
for(i=0;i<=4;i++){
a_x=0;
for(j=0;j<=4;j++) a_x=a_x+a1[i][j]*x[j];
r[i]=b1[i]-a_x;
}
y[0]=r[0];
for(i=2;i<=n;i++){
l_u=0;
for(t=1;t<=i-1;t++)l_u=l_u+l[i-1][t-1]*y[t-1];
y[i-1]=r[i-1]-l_u;
}
d[n-1]=y[n-1]/u[n-1][n-1];
for(i=n-1;i>=1;i--){
u_x=0;
for(t=i+1;t<=n;t++)u_x=u_x+u[i-1][t-1]*d[t-1];
d[i-1]=(y[i-1]-u_x)/u[i-1][i-1];
}
for(i=0;i<=4;i++)x[i]=x[i]+d[i];
}while(MAX(d)/MAX(x)>1e-12);//精度1e-12时中止循环
return 1;
}
void m_Multi(double *a, double *b, double *c, int la, int lb, int lc, int r, int s, int t)
//求r*s阶矩阵A与s*t阶矩阵B的乘积矩阵C
{
int i, j, k;
for (i=0; i<r; i++)
for (j=0; j<t; j++){
*(c+i*lc+j)=0;
for (k=0; k<s; k++)*(c+i*lc+j)+=*(a+i*la+k)*(*(b+k*lb+j));
}
}
void Transpose(double *a, double *ta, int la, int lta, int n, int m)
//求n*m阶矩阵A的转置TA
{
int i, j;
for (i=0; i<m; i++)
for (j=0; j<n; j++) *(ta+i*lta+j)=*(a+j*la+i);
}
double Inverse(double *a, double *b, int la, int lb, int n)
//求n阶方阵A的逆矩阵B
{
int i, j, k;
double temp;
for(i=0; i<n; i++)
for(j=0; j<n; j++)
if (i==j)
*(b+i*lb+j)=1;
else
*(b+i*lb+j)=0;
for (k=0; k<n; k++){
j=k;
for (i=k+1; i<n; i++)
if (fabs(*(a+i*la+k))>fabs(*(a+j*la+k)))j=i;
if (j!=k)
for (i=0; i<n; i++){
temp=*(a+j*la+i);
*(a+j*la+i)=*(a+k*la+i);
*(a+k*la+i)=temp;
temp=*(b+j*lb+i);
*(b+j*lb+i)=*(b+k*lb+i);
*(b+k*lb+i)=temp;
}
if (*(a+k*la+k)==0)
return 0;
if ((temp=*(a+k*la+k))!=1)
for (i=0; i<n; i++){
*(a+k*la+i)/=temp;
*(b+k*lb+i)/=temp;
}
for (i=0; i<n; i++)
if ((*(a+i*la+k)!=0) && (i!=k)){
temp=*(a+i*la+k);
for (j=0; j<n; j++){
*(a+i*la+j)-=temp*(*(a+k*la+j));
*(b+i*lb+j)-=temp*(*(b+k*lb+j));
}
}
}
return 0;
}
double x[5];
double MAX(double a[5])//求最大值
{
int i;
double b[5],max;
for(i=0;i<=4;i++)b[i]=fabs(a[i]);
max=b[0];
for(i=0;i<=4;i++)if(b[i]>max)max=b[i];
return(max);
}
void InLinearEquation()
//求n*m阶矩阵A的转置TA
{
double bb[5],z[5][31][21],dd[5];
int i,j,m;
double jj[5][5],d[5],cc[5],ss[5];
FILE *fp;
for(i=0;i<=30;i++){
for(j=0;j<=20;j++){
z[0][i][j]=6.7;
z[1][i][j]=15.6765-0.270079*(i+j);
z[2][i][j]=5.30249+0.1410822*(i+j);
z[3][i][j]=15;
z[4][i][j]=29.9983-0.2366318*(i+j);
do{
cc[0]=z[0][i][j];
cc[1]=z[1][i][j];
cc[2]=z[2][i][j];
cc[3]=z[3][i][j];
cc[4]=z[4][i][j];
jj[0][0]=-exp(-z[0][i][j]);
jj[1][0]=-2*exp(-2*z[0][i][j]);
jj[2][0]=1+0.008*i;
jj[3][0]=2;
jj[4][0]=1;
jj[0][1]=-2*exp(-2*z[1][i][j]);
jj[1][1]=-exp(-z[1][i][j]);
jj[2][1]=3;
jj[3][1]=1+0.008*j;
jj[4][1]=-2;
jj[0][2]=1;
jj[1][2]=-2;
jj[2][2]=-exp(-z[2][i][j]);
jj[3][2]=1;
jj[4][2]=-3*(1+0.008*i);
jj[0][3]=-2;
jj[1][3]=1;
jj[2][3]=0;
jj[3][3]=exp(-z[3][i][j]);
jj[4][3]=-2*exp(-2*z[3][i][j]);
jj[0][4]=1;
jj[1][4]=-1;
jj[2][4]=-3;
jj[3][4]=-4*exp(-2*z[4][i][j]);
jj[4][4]=-3*exp(-z[4][i][j]);
bb[0]=-(exp(-z[0][i][j])+exp(-2*z[1][i][j])+z[2][i][j]-2*z[3][i][j]+(1+0.008*i)*z[4][i][j]-5.3);
bb[1]=-(exp(-2*z[0][i][j])+exp(-z[1][i][j])-2*z[2][i][j]+(1+0.008*j)*z[3][i][j]-z[4][i][j]+25.6);
bb[2]=-((1+0.008*i)*z[0][i][j]+3*z[1][i][j]+exp(-z[2][i][j])-3*z[4][i][j]+37.8);
bb[3]=-(2*z[0][i][j]+(1+0.008*j)*z[1][i][j]+z[2][i][j]-exp(-z[3][i][j])+2*exp(-z[4][i][j])-31.3);
bb[4]=-(z[0][i][j]-2*z[1][i][j]-3*(1+0.008*i)*z[2][i][j]+exp(-2*z[3][i][j])+3*exp(-z[4][i][j])+42.1);
for(m=0;m<=4;m++){
ss[m]=MAX(jj[m]);
jj[m][m]=jj[m][m]/ss[m];
bb[m]=bb[m]/ss[m];
}
xDooLittle(5,jj,bb);
dd[0]=z[0][i][j]=x[0]+z[0][i][j];
dd[1]=z[1][i][j]=x[1]+z[1][i][j];
dd[2]=z[2][i][j]=x[2]+z[2][i][j];
dd[3]=z[3][i][j]=x[3]+z[3][i][j];
dd[4]=z[4][i][j]=x[4]+z[4][i][j];
d[0]=z[0][i][j]-cc[0];
d[1]=z[1][i][j]-cc[1];
d[2]=z[2][i][j]-cc[2];
d[3]=z[3][i][j]-cc[3];
d[4]=z[4][i][j]-cc[4];
}while(MAX(d)/MAX(cc)>1e-12); //精度1e-12时中止循环
}
}
int k, r, s, p, q, bestr, bests;
double t1[31][31], t2[31][31], t3[31][31], c[9][9],bestc[9][9];
double b[31][9], tb[9][31], g[21][9], tg[9][21];
double temp, acc, best;
if ((fp=fopen("output2.txt", "wb"))==NULL){
printf("Fail to open output file!");
exit(-1);
}
fprintf(fp, "部分数表:%c\n%c\n", 13, k+1, 13);
for (i=0; i<31; i++)
for (j=0; j<9; j++)
b[i][j]=exp(j*log(1+0.008*i));
Transpose(b[0], tb[0], 9, 31, 31, 9);
for (i=0; i<21; i++)
for (j=0; j<9; j++)
g[i][j]=exp(j*log(1+0.008*i));
Transpose(g[0], tg[0], 9, 21, 21, 9);
for (k=0; k<5; k++) {
best=1e-2;
for (r=1; r<=8; r++)
for (s=1; s<=8; s++)
{
m_Multi(tb[0], b[0], t1[0], 31, 9, 31, r+1, 31, r+1);
Inverse(t1[0], c[0], 31, 9, r+1);
m_Multi(c[0], tb[0], t1[0], 9, 31, 31, r+1, r+1, 31);
m_Multi(t1[0], z[k][0], t2[0], 31, 21, 31, r+1, 31, 21);
m_Multi(tg[0], g[0], t1[0], 21, 9, 31, s+1, 21, s+1);
Inverse(t1[0], c[0], 31, 9, s+1);
m_Multi(c[0], tg[0], t3[0], 9, 21, 31, s+1, s+1, 21);
Transpose(t3[0], t1[0], 31, 31, s+1, 21);
m_Multi(t2[0], t1[0], c[0], 31, 31, 9, r+1, 21, s+1);
//求最小二乘法中的拟合系数矩阵C
acc=0;
//求在指定的r和s下,拟合公式的精度
for (i=0; i<31; i++)
for (j=0; j<21; j++){
temp=0;
for (p=0; p<r+1; p++)
for (q=0; q<s+1; q++)
temp+=c[p][q]*b[i][p]*g[j][q];
acc+=(z[k][i][j]-temp)*(z[k][i][j]-temp);
if ((acc>best) && (k!=0)) {
i=30;
j=20;
continue;
}
}
//输出zi=fi(x, y)拟合过程的r, s值和对应精度
if (k==0)fprintf(fp, "%d%c%d%c%.12le%c\n", r, 9, s, 9, acc, 13);
if (acc<best) { //保存最优精度结果
best=acc;
memcpy(bestc, c, sizeof(c));
bestr=r;
bests=s;
}
}
//输出最优拟合方案和对应的精度
fprintf(fp, "Pt(x,y)的系数:"%c\n%c\n", 13, k+1, 13);
fprintf(fp, "%c\nz%d:%c\n", 13, k+1, 13);
fprintf(fp, "k=%d l=%d accuracy: %.12le", bestr, bests, best);
if (best<=1e-4)
fprintf(fp, "<=1e-4");
else
fprintf(fp, ">1e-4 (failed)");
fprintf(fp, "%c\nc(r, s):%c\n", 13, 13);
for (i=0; i<bestr+1; i++){
for (j=0; j<bests+1; j++)
fprintf(fp, "%.12le%c", bestc[i][j], 9);
fprintf(fp, "%c\n", 13);
}
}
printf("Done.\n");
printf("All tasks completed. Open output2.txt to browse the results.\n");
printf("please press Enter to exit!\n");
while (!getch());
fclose(fp);
}
task2输出结果：
部分数表:
1 1 8.837745844330e-001
1 2 6.977120301441e-001
1 3 6.975657242540e-001
1 4 6.999987448592e-001
1 5 1.505611588141e+003
1 6 4.512941773705e+004
1 7 1.872649803994e+005
1 8 8.670652930462e+005
2 1 1.898010526061e-001
2 2 3.280174823932e-003
2 3 3.046761279032e-003
2 4 5.471324216509e-003
2 5 1.504958224150e+003
2 6 4.512987625190e+004
2 7 1.872691715069e+005
2 8 8.670865903482e+005
3 1 1.876627568285e-001
3 2 8.984274352364e-004
3 3 6.274730083732e-004
3 4 3.047559417101e-003
3 5 1.504956026712e+003
3 6 4.512987685427e+004
3 7 1.872691874737e+005
3 8 8.670866404649e+005
4 1 1.873706951595e-001
4 2 5.394108140757e-004
4 3 2.552630812187e-004
4 4 1.683926698549e-003
4 5 1.505543018260e+003
4 6 4.512012150653e+004
4 7 1.872376121115e+005
4 8 8.669618613543e+005
5 1 1.903560150663e+001
5 2 1.884888271124e+001
5 3 1.884848162002e+001
5 4 1.843498381677e+001
5 5 1.769558798928e+003
5 6 4.101866467761e+004
5 7 1.738837249515e+005
5 8 8.140227331282e+005
6 1 1.675253908950e+004
6 2 1.675245017068e+004
6 3 1.675244909018e+004
6 4 1.674980919948e+004
6 5 1.890974466200e+004
6 6 7.448958657218e+003
6 7 1.826432783547e+003
6 8 5.334043019742e+003
7 1 2.333469311487e+003
7 2 2.333299309037e+003
7 3 2.333298145718e+003
7 4 2.329999773526e+003
7 5 5.655699643879e+003
7 6 1.106707317870e+004
7 7 6.809671642811e+004
7 8 3.782409390890e+005
8 1 9.696359203662e+003
8 2 9.696230860322e+003
8 3 9.696229609318e+003
8 4 9.692629298974e+003
8 5 1.280172133745e+004
8 6 2.888465168597e+003
8 7 1.397167732040e+004
8 8 1.063996410661e+005
z1:
k=4 l=3 accuracy: 2.552630812187e-004>1e-4
c(r, s):
-1.475771674723e+005 4.223996105452e+005 -4.029252308762e+005 1.281130833526e+005
5.414906717227e+005 -1.549750705807e+006 1.478313319854e+006 -4.700683934399e+005
-7.440916659950e+005 2.129667465746e+006 -2.031644587878e+006 6.460721025574e+005
4.539439586022e+005 -1.299324668567e+006 1.239638326950e+006 -3.942466781658e+005
-1.037444357402e+005 2.969736562723e+005 -2.833591242344e+005 9.012562268687e+004
z2:
k=4 l=3 accuracy: 2.019681005120e-003>1e-4
c(r, s):
5.023754577229e+005 -1.438698297897e+006 1.373317787138e+006 -4.369770761016e+005
-1.842917169676e+006 5.278217430674e+006 -5.038549373515e+006 1.603296477792e+006
2.532653371180e+006 -7.253724039907e+006 6.924555075671e+006 -2.203548934793e+006
-1.545099847274e+006 4.425390957295e+006 -4.224796867913e+006 1.344518515976e+006
3.530868499781e+005 -1.011343188483e+006 9.655714696860e+005 -3.073129216001e+005
z3:
k=4 l=3 accuracy: 8.616686622126e-004>1e-4
c(r, s):
-3.049601536690e+005 8.732417153219e+005 -8.334163447258e+005 2.651234355719e+005
1.118931421082e+006 -3.204262216118e+006 3.058249273696e+006 -9.729202803344e+005
-1.538009406740e+006 4.404458985904e+006 -4.203838337726e+006 1.337425573074e+006
9.385209624035e+005 -2.687729225892e+006 2.565402509364e+006 -8.162160246987e+005
-2.145266929748e+005 6.143818012870e+005 -5.864528750963e+005 1.866011823656e+005
z4:
k=3 l=3 accuracy: 3.232276740079e-004>1e-4
c(r, s):
-7.981428261743e+003 2.302433528161e+004 -2.222957417791e+004 7.158933779802e+003
2.200807464150e+004 -6.355998062555e+004 6.142910214419e+004 -1.979220689567e+004
-2.016932409158e+004 5.842106931024e+004 -5.652356270043e+004 1.822019003321e+004
6.175464093908e+003 -1.791894243390e+004 1.734385483914e+004 -5.591002131627e+003
z5:
k=4 l=3 accuracy: 2.735948933337e-003>1e-4
c(r, s):
4.417669554680e+005 -1.264908236395e+006 1.207271306427e+006 -3.840995498514e+005
-1.620217650657e+006 4.639831496617e+006 -4.428632496591e+006 1.409072117292e+006
2.226271714202e+006 -6.375467658814e+006 6.085460830414e+006 -1.936343594474e+006
-1.357981813204e+006 3.889006193294e+006 -3.712321301291e+006 1.181322507485e+006
3.102798697371e+005 -8.886322964983e+005 8.483289049685e+005 -2.699773011676e+005
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -