计算实习大作业三
一、算法设计方案：
1.使用牛顿迭代法(简单迭代法不
收敛)，对原题中给出的，，()的
11*21 组分别求出原题中方程组的一组解，于是得到一组和对应的。
2.对于已求出的，使用分片二次代数插 值法对原题中关于的数表进行插值得到。
于是产生了 z=f(x,y)的 11*21 个数值 解。
3.从 k=1 开始逐渐增大 k 的值，并使
用最小二乘法曲面拟合法对 z=f(x,y)进行拟合，得到每次的。当时结束计算，输出拟合
结果。
4. 计 算 的 值 并 输
出结果，以观察逼
近的效果。其中。
二、源代码：
/******************************************
* 计算实习大作业第三题源程序
* 本程序在 Dev-C++ 4.9.9.0 版本下调试通过
* 最后修改：2008.06.11
******************************************/
#define N 6
//矩阵的阶；
#define M 4 //方程组
未知数个数；
#define EPSL 1.0e-12 // 迭代的
精度水平；
#define MAXDIGIT 1.0e-219
#define SIGSUM 1.0e-7
#define L 500 // 迭
代最大次数；
#define K 10 // 最小二乘法拟合时的
次数上限；
#define X_NUM 11
#define Y_NUM 21
#define OUTPUTMODE 1 //输出格式：0--输出至屏幕，1--输出至
文件；
#include <stdio.h>
#include <math.h>
double mat_u[N] = {0, 0.4, 0.8, 1.2, 1.6, 2.0},
mat_t[N] = {0, 0.2, 0.4, 0.6, 0.8, 1.0},
mat_ztu[N][N] = {
{-0.5, -0.34, 0.14, 0.94, 2.06, 3.5},
{-0.42, -0.5, -0.26, 0.3, 1.18, 2.38},
{-0.18, -0.5, -0.5, -0.18, 0.46, 1.42},
{0.22, -0.34, -0.58, -0.5, -0.1, 0.62},
{0.78, -0.02, -0.5, -0.66, -0.5, -0.02},
{1.5, 0.46, -0.26, -0.66, -0.74, -0.5},
},
mat_val_z[X_NUM][Y_NUM] = {0},
init_tuvw[4] = {1,2,1,2},
mat_crs[K][K] = {0};
FILE *p;
int main()//主程序入口
{
int i, j, x, y, kmin, tmp = 0;
double tmpval;
int getval_z(double, double, int, int);
int leasquare();
void result_out(int);
if(OUTPUTMODE)
{
p = fopen("c:\\Result.txt", "w+");
fprintf(p, "计算结果：\n");
fclose(p);
}
for(i = 0; i < X_NUM; i++)
for(j = 0; j < Y_NUM; j++) tmp+= getval_z(0.08 * i,0.5 +
0.05 * j, i, j);
if(!tmp) printf("迭代求解 z=f(x,y) 完毕。\n");
else printf("迭代超过最大次数，计算结果可能不准确！\n");
if(kmin=leasquare())
{
printf("近似表达式已求出！\n");
for(i = 0; i < 8; i++)
for(j = 0; j < 5; j++) tmp+= getval_z(0.1 *
(i+1),0.5 + 0.2 * (j+1), i, j);
if(!tmp)
{
printf("迭代求解 z=f(x*,y*) 完毕。\n");
for(x = 0; x < 8; x++)
{
for(y = 0; y < 5; y++)
{
tmpval = 0;
for(i = 0; i < kmin; i++)
{
for(j = 0; j < kmin; j++)
{
tmpval= tmpval + mat_crs[i][j] * pow(0.1
* (x+1), i) * pow(0.5 + 0.2 * (y+1), j);
}
}
if(OUTPUTMODE)
{
p = fopen("c:\\Result.txt", "at+");
fprintf(p, "x*[%d]=%f, y*[%d]=%f: f(x*[%d],
y*[%d])=%.12le, p(x*[%d], y*[%d])=%.12le\n", x+1, 0.1 * (x+1),
y+1, 0.5 + 0.2 * (y+1), x+1, y+1, mat_val_z[x][y], x+1, y+1,
tmpval);
fclose(p);
}
else
printf("x*[%d]=%f, y*[%d]=%f: f(x*[%d], y*[%d])=
%.12le, p(x*[%d], y*[%d])=%.12le\n", x+1, 0.1 * (x+1), y+1, 0.5 +
0.2 * (y+1), x+1, y+1, mat_val_z[x][y], x+1, y+1, tmpval);
}
}
}
printf("结果见 C:\\Result.txt！");
getchar();
}
else
{
printf("近似表达式未求出，请增大 K 值后重试！\n");
getchar();
}
}
int getval_z(double x, double y, int i, int j)//使用牛顿迭代法迭代求
解方程组 f(x,y,t,u,v,w)。
{
double vec_tuvw[M], vec_delta[M], fdao[M][M], f[M], mo_k,
mo_delta_k, shang, t, u;
int n = 0, k, flag_max, flag_stat;
void solve(double fdao[M][M], double vec_delta[M], double f[M]);
double interpolation(double, double);
flag_stat = 1;
for(k = 0; k < M; k++) vec_tuvw[k] = init_tuvw[k];
do
{
f[0] = -1.0 *(0.5 * cos(vec_tuvw[0]) + vec_tuvw[1] +
vec_tuvw[2] + vec_tuvw[3] - x - 2.67);
f[1] = -1.0 *(vec_tuvw[0] + 0.5 * sin(vec_tuvw[1]) +
vec_tuvw[2] + vec_tuvw[3] - y - 1.07);
f[2] = -1.0 *(0.5 * vec_tuvw[0] + vec_tuvw[1] +
cos(vec_tuvw[2]) + vec_tuvw[3] - x - 3.74);
f[3] = -1.0 *(vec_tuvw[0] + 0.5 * vec_tuvw[1] +
vec_tuvw[2] + sin(vec_tuvw[3]) - y - 0.79);
fdao[0][0] = -0.5 * sin(vec_tuvw[0]);
fdao[0][1] = 1;
fdao[0][2] = 1;
fdao[0][3] = 1;
fdao[1][0] = 1;
fdao[1][1] = 0.5 * cos(vec_tuvw[1]);
fdao[1][2] = 1;
fdao[1][3] = 1;
fdao[2][0] = 0.5;
fdao[2][1] = 1;
fdao[2][2] = -1 * sin(vec_tuvw[2]);
fdao[2][3] = 1;
fdao[3][0] = 1;
fdao[3][1] = 0.5;
fdao[3][2] = 1;
fdao[3][3] = cos(vec_tuvw[3]);
solve(fdao, vec_delta, f);
flag_max = 0;
for(k = 1; k < M; k++) if(fabs(vec_delta[k]) >
fabs(vec_delta[flag_max])) flag_max = k;
mo_delta_k = vec_delta[flag_max];
flag_max = 0;
for(k = 1; k < M; k++) if(fabs(vec_tuvw[k]) >
fabs(vec_tuvw[flag_max])) flag_max = k;
mo_k = vec_tuvw[flag_max];
shang = fabs(mo_delta_k / mo_k);
if(shang < EPSL)
{
t = vec_tuvw[0] + vec_delta[0], u = vec_tuvw[1] +
vec_delta[1];
flag_stat = 0;
break;
}
else
for(k = 0; k < M; k++) vec_tuvw[k]+= vec_delta[k];
}
while(n++ <= L);
if(!flag_stat)
{
if(OUTPUTMODE)
{
p = fopen("c:\\Result.txt", "at+");
fprintf(p, "x%d=%f, y%d=%f: f(x%d, y%d)=%.12le\n", i,
0.08*i, j, 0.5+0.05*j, i, j, mat_val_z[i][j] =
interpolation(u,t));
fclose(p);
}
else
printf("x%d=%f, y%d=%f: f(x%d, y%d)=%.12le\n", i, 0.08*i,
j, 0.5+0.05*j, i, j, mat_val_z[i][j] = interpolation(u,t));
}
return flag_stat;
}
void solve(double fdao[M][M], double vec_delta[M], double f[M])//
牛顿消元法解线性方程组；
{
int i, j, k, ik;
double tmp, mik;
for(k = 0; k < M-1; k++)
{
ik = k;
for(i = k; i < M; i++) if(fabs(fdao[ik][k]) <
fabs(fdao[i][k])) ik = i;
for(j = k; j < M; j++)
- 1
- 2
- 3
前往页