对于龙格函数 f ( x ) = 1 ( 1 + 25 x 2 ) f(x)=\frac{1}{(1+25x^2)} f(x)=(1+25x2)1在区间 [ − 1 , 1 ] [-1,1] [−1,1]上取 x i = − 1 + 0.2 i ( i = 0 , 1 , . . , 10 ) x_i=-1+0.2i(i=0,1,..,10) xi=−1+0.2i(i=0,1,..,10),试求3次曲线拟合,试画出拟合曲线并打印出方程。
一、C语言代码
#include<stdio.h>
#include<math.h>
#define N 11
#define M 4
int main()
{
float x[N] = {
0 }, y[N] = {
0 }, Matrix[M][M], B[M] = {
0 },a[M];
float sum,temp_1,temp_2;
int i,j,k;
printf("自变量x取值(11个):\n");
for (i = 0; i < 11; i++)
{
x[i] = -1 + 0.2*i;
printf("%.6lf ", x[i]);
}
printf("\n\n");
printf("因变量y取值(11个):\n");
for (i = 0; i < 11; i++)
{
y[i] = 1 / (1 + 25 * pow(x[i], 2));
printf("%.6lf ", y[i]);
}
printf("\n\n");
for (i = 0; i < M; i++)
{
for (j = i; j < M; j++)
{
sum = 0.0;
for (k = 0; k < 11; k++)
{
sum += pow(x[k],i) * pow(x[k],j);
}
Matrix[i][j] = Matrix[j][i] = sum;
}
}
printf("输入系数矩阵Matrix:\n");
for (i = 0; i < M; i++)
{
for (j = 0; j < M; j++)
{
printf("%.6lf\t", Matrix[i][j]);
}
printf("\n");
}
printf("\n");
for (i = 0; i < M; i++)
{
sum = 0.0;
for (k = 0; k < 11; k++)
{
sum += pow(x[k], i) * y[k];
}
B[i] = sum;
}
printf("右端的输入项B:\n");
for (i = 0; i < M; i++)
{
printf("%.6lf\n", B[i]);
}
printf("\n");
for (k = 0; k < M - 1; k++)
{
if (!Matrix[k][k])
{
return -1;
}
for (i = k + 1; i < M; i++)
{
temp_1 = Matrix[i][k] / Matrix[k][k];
for (j = k; j < M; j++)
{
Matrix[i][j] -= temp_1 * Matrix[k][j];
}
B[i] -= temp_1 * B[k];
}
}
printf("消元后的矩阵A:\n");
for (i = 0; i < M; i++)
{
for (j = 0; j < M; j++)
{
printf("%.6lf\t", Matrix[i][j]);
}
printf("\n");
}
printf("\n");
printf("消元后的右端B:\n");
for (i = 0; i < M; i++)
{
printf("%.6lf\t", B[i]);
printf("\n");
}
a[M - 1] = B[M - 1] / Matrix[M - 1][M - 1];
for (k = M - 2; k >= 0; k--)
{
temp_2 = B[k];
for (j = k + 1; j < M; j++)
{
temp_2 -= Matrix[k][j] * a[j];
}
a[k] = temp_2 / Matrix[k][k];
printf("\n");
}
printf("结果a=:\n");
for (i = 0; i < M; i++)
{
printf("a%d = %.6lf\n",i,a[i]);
}
printf("\n");
printf("三次拟合多项式方程为:\n");
printf("P(x)=%.6lf+%.6lfx%.6lfx^2%.6lfx^3",a[0], a[1], a[2], a[3]);
printf("\n\n");
printf("完善后三次拟合多项式方程为:\n");
printf("P(x)=%.6lf%.6lfx^2:\n\n", a[0], a[2]);
return 0;
}