数值分析实习作业三

数值分析大作业第三题学号：姓名：学院：授课教师：吕淑娟2016.11.21题目：关于x, y, t, u, v, w的下列方程组（A.3） (A.3)以及关于z, t, u的下列二维数表z ut00.40.81.21.62.00-0.5-0.340.140.942.063.50.2-0.42-0.5-0.260.31.182.380.4-0.18-0.5-0.5-0.180.461.420.60.22-0.34-0.58-0.5-0.10.620.80.78-0.02-0.5-0.66-0.5-0.021.01.50.46-0.26-0.66-0.74-0.5确定了一个二元函数z = f(x, y)。1、试用数值方法求出f(x, y)在区域 上的一个近似表达式要求以最小的k值达到以下的精度其中，。2、计算（i = 1, 2, ，8；j = 1, 2,，5）的值，以观察逼近的效果，其中，=0.1i, =0.5+0.2j。说明：1、用迭代方法求解非线性方程组时，要求近似解向量满足2、作二元插值时，要使用分片二次代数插值。3、要由程序自动确定最小的k值。4、打印以下内容：（1）算法的设计方案；（2）全部源程序； （3）数表： （i = 0,1,2,10；j = 0,1,2,20）；（4）选择过程的值；（5）达到精度要求时的值以及中的系数（r = 0,1,k；s = 0,1,k）；（6）数表：（i = 1, 2, ,8；j = 1, 2,5）；（7）讨论：计算过程中遇到的现象。5、采用f型输出的准确值，其余实型数采用e型输出并且至少显示12位有效数字。一、 算法设计方案1、 方程组(A.3)由四个方程组成，又为六元方程组，则给x,y赋一组值，则能解出一组（t,u,v,w）。依次代入xi=0.08i，yj=0.5+0.05j (其中i=0,1,10; j=0,1,20)，利用NEWTON迭代法解得对应的(tij,uij)。NEWTON迭代法具体算法见教材p9091。2、 利用二维数表进行分片二次代数插值，得到函数z=z(t,u)。对于区域D内任意x,y，代入方程组（A.3）解得(t,u)，再代入函数z=z(t,u)，此从x,y到z的映射即为二元函数 z=f(x,y)。分片二次代数插值具体算法见教材p101。3、 将1中得到的(tij,uij)代入函数z=z(t,u)，得到zij=f(xi,yj)。利用最小二乘法进行拟合，得到，其中使k依次取1,2,直到使拟合函数满足。则所取k为满足条件的最小值。最小二乘法拟合具体算法见教材p142。4、 同1、2步，代入=0.1i, =0.5+0.2j，（i = 1, 2, ，8；j= 1, 2,，5），得到；将代入 ，得到。比较两函数值，观察逼近的效果。二、 结果1、数表： （i = 0,1,2,10；j = 0,1,2,20）2、选择过程的值3、最终k应为5，此时crs矩阵如下4、数表：（i = 1, 2, ,8；j = 1, 2,5）三、 遇到的问题从插值得到的f(x,y)函数开始（包括其在内），所得到计算结果与课本答案存在1e-9左右的差异。经验证，并不是迭代初值的选择所导致的。四、 源程序#include stdafx.h#include stdio.h#include math.h#include stdlib.h#define N 4#define L 11#define M 21int _tmain(int argc, _TCHAR* argv)void ntsl(double x, double y, double aN);/Newton迭代解方程组函数void gsslm(double aMM, double bM, double xM, int n);void gssln(double aNN, double bN, double xN);/列主元gauss迭代void fcfz(double aN, double fN, double x, double y);/方程组赋值函数void fcds(double aN, double ffNN);/方程导数矩阵函数double amul(double aN, double bN);/向量内积函数double fpcz(double t03, double u03, double z033, double t, double u);/二次分片插值函数void atra(double aMM, int n, int m, double atMM);/矩阵转置void aamu(double aMM, int n, int m, int s, double bMM, double cMM);/矩阵相乘void aqni(double aMM, int n, double aniMM);/矩阵求逆double pxy(double x, double y, double cMM, int k);/求p(x,y)void zxec(int l, int m, double uLM, double xL, double yM, double cMM);/最小二乘double ztu(double t, double u, double z0066);/求z(t,u)double tuvwN = 1,1,1,1;/方程组迭代赋初值double x;double y;double zLM;int k0;double cMM = 0 ;double x1L;double y1M;double z0066 = -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 ;for (int i = 0; i L; i+)for (int j = 0; j 7; j+)if (i = 0)if (j = 0)printf(x-y );printf( %1.2lf , 0.5 + 0.05*j);if (j = 6)printf(n);for (int j = 0; j 7; j+)x = 0.08*i;y = 0.5 + 0.05*j;ntsl(x, y, tuvw);zij = ztu(tuvw0, tuvw1,z00);if (j = 0)printf(%1.2lf , 0.08*i);printf(% 1.11e , zij);if (j = 6)printf(n);if (i = 10)printf(n);for (int i = 0; i L; i+)for (int j = 7; j 14; j+)if (i = 0)if (j = 7)printf(x-y );printf( %1.2lf , 0.5 + 0.05*j);if (j = 13)printf(n);for (int j = 7; j 14; j+)x = 0.08*i;y = 0.5 + 0.05*j;ntsl(x, y, tuvw);zij = ztu(tuvw0, tuvw1, z00);if (j = 7)printf(%1.2lf , 0.08*i);printf(% 1.11e , zij);if (j = 13)printf(n);if (i = 10)printf(n);for (int i = 0; i L; i+)for (int j = 14; j 21; j+)if (i = 0)if (j = 14)printf(x-y );printf( %1.2lf , 0.5 + 0.05*j);if (j = 20)printf(n);for (int j = 14; j 21; j+)x = 0.08*i;y = 0.5 + 0.05*j;ntsl(x, y, tuvw);zij = ztu(tuvw0, tuvw1, z00);if (j = 14)printf(%1.2lf , 0.08*i);printf(% 1.11e , zij);if (j = 20)printf(n);if (i = 10)printf(n);for (int i = 0; i L; i+)x1i = 0.08*i;for (int i = 0; i M; i+)y1i = 0.5 + 0.05*i;/最小二乘法拟合+循环确定满足条件k值for (int k = 0; k L - 1; k+)double sum1LM = 0 ;double sum2 = 0;zxec(k, k, z, x1, y1, c);for (int i = 0; i L; i+)for (int j = 0; j M; j+)for (int r = 0; r k + 1; r+)for (int s = 0; s k + 1; s+)if (r=0) & (i=0)sum1ij = sum1ij + crs * 1 * pow(y1j, s);elsesum1ij = sum1ij + crs * pow(x1i, r)*pow(y1j, s);sum2 = sum2 + pow(sum1ij-zij),2);printf(k=%d , k);printf(=%1.11en,sum2);if (sum2 (1e-7)printf(k=%dn, k);k0 = k;break;if (k = L - 2)printf(拟合失败n);for (int i = 0; i 6; i+)for (int j = 0; j 6; j+)printf(% .11e , cij);if (j = 5)printf(n);for (int i = 1; i 9; i+)/求f(x*,y*)if (i = 1)printf(x-y-f );for (int j = 1; j 6; j+)if (i = 1)printf( %.2lf , 0.5 + 0.2*j);if (j = 5)printf(n);printf(%.2lf , 0.1*i);for (int j = 1; j 6; j+)x = 0.1*i;y = 0.5 + 0.2*j;ntsl(x, y, tuvw);printf(% .11e , ztu(tuvw0, tuvw1, z00);if (j = 5)printf(n);printf(n);/求p(x*,y*)for (int i = 1; i 9; i+)if (i = 1)printf(x-y-p );for (int j = 1; j 6; j+)if (i = 1)printf( %.2lf , 0.5 + 0.2*j);if (j = 5)printf(n);printf(%.2lf , 0.1*i);for (int j = 1; j 6; j+)x = 0.1*i;y = 0.5 + 0.2*j;printf(% .11e , pxy(x, y, c, k0);if (j = 5)printf(n);return 0; void gssln(double aNN, double bN, double xN)double aiN;double bi;for (int i = 0; i N; i+)xi = 0;for (int k = 0; k N; k+)double max = 0;int s = 0;double mN = 0 ;for (int i = k; i N; i+)if (max fabs(aik)max = fabs(aik);s = i;for (int i = 0; i N; i+)aii = aki;aki = asi;asi = aii;bi = bk;bk = bs;bs = bi;for (int i = k + 1; i N; i+)mi = aik / akk;for (int j = k + 1; j -1; k-)for (int j = k + 1; j N; j+)xk = xk + xj * akj;xk = bk / akk - xk/akk;void fcfz(double aN, double fN, double x, double y)f0 = cos(a0) / 2 + a1 + a2 + a3 - 2.67 - x;f1 = a0 + sin(a1) / 2 + a2 + a3 - 1.07 - y;f2 = a0 / 2 + a1 + cos(a2) + a3 - 3.74 - x;f3 = a0 + a1 / 2 + a2 + sin(a3) - 0.79 - y;void fcds(double aN, double ffNN)for (int i = 0; i N; i+)for (int j = 0; j N; j+)ffij = 1;ff00 = 0 - sin(a0) / 2;ff11 = cos(a1) / 2;ff20 = 0.5;ff22 = 0 - sin(a2);ff31 = 0.5;ff33 = cos(a3);void ntsl(double x, double y, double aN)double fN = 0 ;double ffNN;double detaN = 0 ;int q = 0;for (int s = 0; s 10000; s+)fcfz(a, f, x, y);fcds(a, ff);double f1N;for (int i = 0; i N; i+)f1i = 0 - fi;gssln(ff, f1, deta);/求解线性方程组for (int i = 0; i N; i+)ai = ai + detai;for (int i = 0; i 1e-12)break;if (i = 3) q = 1;if (q = 1)break;if (s = 9999)printf(方程迭代失败n);system(pause);double amul(double aN, double bN)double r=0;for (int i = 0; i N; i+)r = r + ai * bi;return r;double fpcz(double t03,double u03, double z033,double t,double u)double lt3 = 1,1,1;double lu3 = 1,1,1;double fi33;double z=0;for (int k = 0; k 3; k+)for (int i = 0; i 3; i+)if (i != k)ltk = ltk * (t - t0i) / (t0k - t0i);luk = luk * (u - u0i) / (u0k - u0i);for (int i=0; i 3; i+)for (int j=0; j 3; j+)fiij = lti * luj;z = z + fiij * z0ij;return z;double ztu(double t, double u,double z0066)double z;int s;int q;double t03;double u03;double z033;if (t 0.9)s = 4;elses = (t + 0.1) / 0.2;for (int r = 0; r3; r+)t0r = 0.2*(s + r - 1);if (u 1.8)q = 4;elseq = (u + 0.2) / 0.4;for (int r = 0; r3; r+)u0r = 0.4*(q + r - 1);for (int r = 0; r 3; r+)for (int k = 0; k 3; k+)z0rk = z00s + r - 1q + k - 1;z = fpcz(t0, u0, z0, t, u);return z;void aamu(double aMM, int n, int m, int s, double bMM, double cMM)for (int i = 0; i M; i+)for (int j = 0; j M; j+)cij = 0;for (int i = 0; i n; i+)for (int j = 0; j s; j+)for (int k = 0; k m; k+)cij = cij + aik * bkj;void atra(double aMM, int n, int m, double atMM)for (int i = 0; i n; i+)for (int j = 0; j m; j+)atji = aij;void gsslm(double aMM, double bM, double xM, int n)double aiM;double bi;for (int i = 0; i M;i+)xi =0;for (int k = 0; k n; k+)double max = 0;int s = 0;double mM = 0 ;for (int i = k; i n; i+)if (max fabs(aik)max = fabs(aik);s = i;for (int i = 0; i n; i+)aii = aki;aki = asi;asi = aii;bi = bk;bk = bs;bs = bi;for (int i = k + 1; i n; i+)mi = aik / akk;for (int j = k + 1; j -1; k-)for (int j = k + 1; j n; j+)xk = xk + xj * akj;xk = bk / akk - xk/akk;void aqni(double aMM, int n, double aniMM)double IMM = 0 ;double anitMM = 0 ;for (int i = 0; i M; i+)Iii = 1;for (int i = 0; i n; i+)double a1MM;for (int i = 0; i n; i+)for (int j = 0; j n; j+)a1ij = aij;gsslm(a1, Ii, aniti, n);atra(anit, n, n, ani);void zxec(int l, int m, double uLM, double xL,double yM,double cMM)double bMM = 0;double gMM = 0;for (int i = 0; i L; i+)for (int j = 0; j l+1; j+)if (i = 0) & (j = 0)bij = 1;else bij = pow(xi, j);for (int i = 0; i M; i+)for (int j = 0; j m + 1; j+)gij = pow(yi, j);double btMM = 0;double gtMM = 0;atra(b, L, l + 1, bt);atra(g, M, m + 1, gt);double btbMM = 0 ;double gtgMM = 0 ;aamu(bt, l + 1,L,l+1, b,btb);aamu(gt,m + 1,M, m+1,g, gtg);aqni(btb, l + 1, b);aqni(gtg, m + 1, gt);aamu(b, l + 1, l + 1, L, bt, btb);aamu(g, M, m + 1, m + 1, gt, gtg);aamu(btb, l + 1, L, M, u, b);aamu(b, l + 1, M, m + 1, gtg, c);double pxy(double x, double y, double cMM, int k)double p=0;for (int r = 0; r k + 1; r+)for (int s = 0; s k + 1; s+)#include stdafx.h#include stdio.h#include math.h#include stdlib.h#define N 4#define L 11#define M 21int _tmain(int argc, _TCHAR* argv)void ntsl(double x, double y, double aN);/Newton迭代解方程组函数void gsslm(double aMM, double bM, double xM, int n);void gssln(double aNN, double bN, double xN);/列主元gauss迭代void fcfz(double aN, double fN, double x, double y);/方程组赋值函数void fcds(double aN, double ffNN);/方程导数矩阵函数double amul(double aN, double bN);/向量内积函数double fpcz(double t03, double u03, double z033, double t, double u);/二次分片插值函数void atra(double aMM, int n, int m, double atMM);/矩阵转置void aamu(double aMM, int n, int m, int s, double bMM, double cMM);/矩阵相乘void aqni(double aMM, int n, double aniMM);/矩阵求逆double pxy(double x, double y, double cMM, int k);/求p(x,y)void zxec(int l, int m, double uLM, double xL, double yM, double cMM);/最小二乘double ztu(double t, double u, double z0066);/求z(t,u)double tuvwN = 1,1,1,1;/方程组迭代赋初值double x;double y;double zLM;int k0;double cMM = 0 ;double x1L;double y1M;double z0066 = -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 ;for (int i = 0; i L; i+)for (int j = 0; j 7; j+)if (i = 0)if (j = 0)printf(x-y );printf( %1.2lf , 0.5 + 0.05*j);if (j = 6)printf(n);for (int j = 0; j 7; j+)x = 0.08*i;y = 0.5 + 0.05*j;ntsl(x, y, tuvw);zij = ztu(tuvw0, tuvw1,z00);if (j = 0)printf(%1.2lf , 0.08*i);printf(% 1.11e , zij);if (j = 6)printf(n);if (i = 10)printf(n);for (int i = 0; i L; i+)for (int j = 7; j 14; j+)if (i = 0)if (j = 7)printf(x-y );printf( %1.2lf , 0.5 + 0.05*j);if (j = 13)printf(n);for (int j = 7; j 14; j+)x = 0.08*i;y = 0.5 + 0.05*j;ntsl(x, y, tuvw);zij = ztu(tuvw0, tuvw1, z00);if (j = 7)printf(%1.2lf , 0.08*i);printf(% 1.11e , zij);if (j = 13)printf(n);if (i = 10)printf(n);for (int i = 0; i L; i+)for (int j = 14; j 21; j+)if (i = 0)if (j = 14)printf(x-y );printf( %1.2lf , 0.5 + 0.05*j);if (j = 20)printf(n);for (int j = 14; j 21; j+)x = 0.08*i;y = 0.5 + 0.05*j;ntsl(x, y, tuvw);zij = ztu(tuvw0, tuvw1, z00);if (j = 14)printf(%1.2lf , 0.08*i);printf(% 1.11e , zij);if (j = 20)printf(n);if (i = 10)printf(n);for (int i = 0; i L; i+)x1i = 0.08*i;for (int i = 0; i M; i+)y1i = 0.5 + 0.05*i;/最小二乘法拟合+循环确定满足条件k值for (int k = 0; k L - 1; k+)double sum1LM = 0 ;double sum2 = 0;zxec(k, k, z, x1, y1, c);for (int i = 0; i L; i+)for (int j = 0; j M; j+)for (int r = 0; r k + 1; r+)for (int s = 0; s k + 1; s+)if (r=0) & (i=0)sum1ij = sum1ij + crs * 1 * pow(y1j, s);elsesum1ij = sum1ij + crs * pow(x1i, r)*pow(y1j, s);sum2 = sum2 + pow(sum1ij-zij),2);printf(k=%d , k);printf(=%1.11en,sum2);if (sum2 (1e-7)printf(k=%dn, k);k0 = k;break;if (k = L - 2)printf(拟合失败n);for (int i = 0; i 6; i+)for (int j = 0; j 6; j+)printf(% .11e , cij);if (j = 5)printf(n);for (int i = 1; i 9; i+)/求f(x*,y*)if (i = 1)printf(x-y-f );for (int j = 1; j 6; j+)if (i = 1)printf( %.2lf , 0.5 + 0.2*j);if (j = 5)printf(n);printf(%.2lf , 0.1*i);for (int j = 1; j 6; j+)x = 0.1*i;y = 0.5 + 0.2*j;ntsl(x, y, tuvw);printf(% .11e , ztu(tuvw0, tuvw1, z00);if (j = 5)printf(n);printf(n);/求p(x*,y*)for (int i = 1; i 9; i+)if (i = 1)printf(x-y-p );for (int j = 1; j 6; j+)if (i = 1)printf( %.2lf , 0.5 + 0.2*j);if (j = 5)printf(n);printf(%.2lf , 0.1*i);for (int j = 1; j 6; j+)x = 0.1*i;y = 0.5 + 0.2*j;printf(% .11e , pxy(x, y, c, k0);if (j = 5)printf(n);return 0; void gssln(double aNN, double bN, double xN)double aiN;double bi;for (int i = 0; i N; i+)xi = 0;for (int k = 0; k N; k+)double max = 0;int s = 0;double mN