数值分析课程作业

数值分析课程作业

（1）解答如下：

不稳定的数值计算公式往往会出现“差之毫厘，失之千里”的错误结果。因此，在计算过程中要选用稳定的计算公式。

在本题的迭代计算中，由I0往后递推的求解I n会使误差快速增大，而先近似计算I

然后回代计算I n-1…I0会使误差逐渐减小，因此使用此计算方法较为稳定。

n

用C语言实现的代码如下：

#include

/* 计算 k的阶乘 */

int k_factorial(int k)

{

if (0 == k)

return 1;

else

return k * k_factorial(k-1);

}

/* 计算 1/e，k控制计算精度 */

double get_e_1(int k)

{

if (0 == k)

return 1;

else

return 1.0*((k%2) ? -1 : 1)/k_factorial(k) + get_e_1(k-1);

}

/*

* 该函数用于计算公式的近似值。对于给定的n值，首先计算(n+k)值对应的

* 公式的近似值，然后回代计算n值对应的近似值。

*/

double get_value(int n, int k)

{

double res;

res = (1+get_e_1(7)) / ((n+k)+1) / 2;

while (k > 0) {

res = (1-res) / (n+k);

k--;

}

return res;

}

int main(int argc, char* argv[]) { int n;

printf("请输入一个n 值："); scanf("%d", &n); if (n < 0)

return 0;

printf("对于给定的n=%d ，公式的近似值为：%lf\n", n, get_value(n, 9)); return 0;

}

（2） 解答如下：

考虑拟合函数：)()()()(1100x q a x q a x q a x n n +++= ϕ，将数据表

⎪⎪

⎩⎪⎪⎨

⎧=++++=++++=++++m m n n m m m n n n n y x q a x q a x q a x q a y x q a x q a x q a x q a y x q a x q a x q a x q a )()()()()()()()()()()()(2211002222221120011122111100 （m > n ）

其系数矩阵为

⎥⎥⎥

⎥

⎦⎤⎢⎢⎢⎢⎣⎡)()()()()()()()()()()()(21022221201121110m n m m m

n n x q x q x q x q x q x q x q x q x q x q x q x q

由于多项式q 0(x )，q 1(x )，q 2(x )，……，q n (x )在点集{x 1，x 2，……，x m }上的正交，所以超定方程组的系数矩阵中不同列的列向量是相互正交的向量组。于是用这一矩阵的转置矩阵去左乘超定方程组左、右两端得正规方程组

⎪⎪⎩⎪⎪⎨⎧===),(),(),(),(),(),(11110000y q a q q y q a q q y q a q q n n n n => ⎪⎪⎩⎪⎪⎨

⎧===),/(),(),/(),(),/(),(11110000n n n n q q y q a q q y q a q q y q a

其中，∑==m

i i k k k x q q q 12

)

(),(，

∑==m

i i

i k k y x q y q 1

)(),(。因为正规方程组中每一个方程都

是一元一次方程可以直接写出原超方程组的最小二乘解，所以拟合函数为

)

(),(),()(),()

,()(),(),()(11110000x q q q y q x q q q y q x q q q y q x n n n n +++=

ϕ

这一结果与用次多项式拟合所得结果在理论是完全一样的，只是形式上不同、算法

实现上避免了解病态方程组。

用C 语言实现的代码如下：

#include #include #include #include

const double EPS=1E-10; // 运算精度 /**

* 读入节点数组x 、函数值数组y 、权值数组ω及节点数N */

void readFile(double *&x,double *&y,double *&omiga,int &N) {

FILE *fp;

fp=fopen(".\\ZXECF.txt","r"); if (fp==NULL) {

printf("指定位置的文件不存在，请检查！！"); getch(); exit(0); }

for (N=1; ;N++) { // 计算节点个数N fscanf(fp,"%*lf%*lf%*lf"); if (feof(fp)) break; }

x=new double[N]; y=new double[N]; omiga=new double[N];

rewind(fp);

for (int i=0;i

fscanf(fp,"%lf%lf%lf",&x[i],&y[i],&omiga[i]);

fclose(fp); } /**