数值分析上机实验报告之样条插值.docx

《数值分析上机实验报告之样条插值.docx》由会员分享，可在线阅读，更多相关《数值分析上机实验报告之样条插值.docx（24页珍藏版）》请在沃文网上搜索。

1、数值分析上机实验报告之样条插值1.三次样条插值（初值条件1）：P52.9、给定函数y=fx的函数表和边界条件s75=0，s80=0，求三次样条插值函数s(x)，并求f(78.3)的近似值。函数表x757677787980y=fx2.7682.8332.9032.9793.0623.153源代码：yangtiao.cpp#include#includevoid main()int choice = 0;int n = 2;double xx,*x, *y,*a,*b,*a1,*b1,*h,*m;cout请输入插值节点个数n：n;x = new doublen; y = new doublen;a

2、 = new doublen; b = new doublen;a1 = new doublen; b1 = new doublen;h = new doublen-1; m = new doublen+1;cout请输入n个插值的节点(xi，yi)：endl;for (int i = 0 ; i xiyi;for (int j = 0 ; j n-1 ; j+)hj = xj+1 - xj;cout请输入待估点xx：xx;cout请选择边界条件：choice;switch(choice)case 1:double temp1,temp2;a0 = 0;an-1 = 1;cout请输入边界条件

3、的两个一阶微商值s(x1)与s(xn)：temp1temp2;b0 = 2*temp1;bn-1 = 2*temp2;break;case 2:a0 = 1;an-1 = 0;b0 = 3/h0*(y1-y0);bn-1 = 3/hn-2*(yn-1-yn-2);break;for (int k = 1 ; kn-1 ; k+)ak = hk-1/(hk-1+hk);bk = 3*( (1-ak)/hk-1*(yk-yk-1) + ak/hk*(yk+1-yk) );a10 = -a0/2;b10 = b0/2;for (int l = 1; l=0; j-)mj = a1j * mj+1 +

4、 b1j;/判别xx所在区间并输出结果coutn插值结果为：;for(k = 0 ;k n-1 ;k+)if (xk xx )double output = 0;output = (1+2*(xx-xk)/(xk+1-xk)* pow(xx-xk+1)/(xk-xk+1),2) *yk +(1+2*(xx-xk+1)/(xk-xk+1)* pow(xx-xk)/(xk+1-xk),2) *yk+1 +(xx-xk)* pow(xx-xk+1)/(xk-xk+1),2) *mk+(xx-xk+1)* pow(xx-xk)/(xk+1-xk),2) *mk+1;cout output endl;b

5、reak;delete x;delete y; delete a; delete b; delete a1;delete b1; delete h; delete m;运行结果截图：2.三次样条插值（初值条件2）：P52.10、给定函数y=fx的函数表和边界条件s0.25=1，s0.53=0.6868，求三次样条插值函数s(x)，并求f(0.35)的近似值。函数表x0.250.30.390.450.53y=fx0.50.54770.62450.67080.728源代码：yangtiao.cpp(同上)运行结果截图：3.自动选取步长梯形法：P977、使用自动选取步长梯形法计算积分0121+t2d

6、t的近似值。（给定=0.01）源代码：SelfSelLength.cpp#include#includedouble fun(double a)return 2/( 1+a*a );double SelfSelLength(double R_a,double R_b,double e)double h = (R_b-R_a)/2;double R1 =(fun(R_a)+fun(R_b) * h;int n = 1;double R0;double S;double E;do /每当误差值不符合要求时，计算下一个result值R0 = R1;S = 0;for (int k =1 ;k 3*e

7、);return R1;void main()double a,b,e;cout 请依次输入待求积分函数的下界a、上界b 及 精度要求e: a b e;cout 自动选取步长梯形法可求得积分: SelfSelLength (a,b,e)endl;运行结果截图：4.Romberg求积法：P978、使用Romberg求积法计算积分19xdx的近似值。（给定=0.01，且取n=1）源代码：Romberg.cpp#include#includestatic double Tri128128;double fun(double a)return sqrt(a);double Romberg(double

8、 R_a,double R_b,double e)Tri00 = (R_b-R_a)/2*(fun(R_a)+fun(R_b);int k = 0;double E;do /每当误差值不符合要求时，计算下一行的Tri值k+;double temp = 0;/计算T0k的数值for (int i = 1 ; i=pow(2,k-1) ;i+)temp = temp + fun(R_a + (2*i-1)* (R_b - R_a)/pow(2,k);Tri0k = 0.5 * (Tri0k-1 + (R_b-R_a)/pow(2,k-1) * temp);for (int m = 1 ; m e)

9、;return Trik0;void main()double a,b,e;cout 请依次输入待求积分函数的下界a、上界b 及 精度要求e: a b e;cout 按Romberg求积法可求得积分: Romberg (a,b,e)endl;运行结果截图：5.列主元高斯消去法：测试矩阵为：2-100-12-100-12-100-12x1x2x3x4 = 1010源代码：Guess_Elimination.cpp/列主元高斯消去法#include#include#define N 4 /矩阵的维数，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,

10、2,-1,0,0,-1,2;/系数矩阵static double BN=1,0,1,0; /右端项static double XN; int i,j,k; /计数器void main()for(k = 0; k N-1 ;k+)/选取最大主元int index = k;for(i = k; i N ;i+)if(fabs(Aindexk) fabs(Aik)index = i;/交换行double temp;for( i = k ; i N ;i+ ) temp = Aindexi; Aindexi = Aki; Aki = temp;temp = Bindex;Bindex = Bk;Bk

11、= temp;for(i = k+1; iN; i+)double T = Aik/Akk;Bi = Bi - T * Bk;for ( j = k+1 ; j =0 ; i-)double Temp = 0;for (int j = i+1; jN ;j+)Temp = Temp + Aij * Xj;Xi = (Bi - Temp) /Aii;cout 线性方程组的解（X1，X2，X3.Xn）为：endl;for( i = 0; i N ;i+)cout Xi ;运行结果截图：6.列主元LU分解法：测试矩阵为：2-100-12-100-12-100-12x1x2x3x4 = 1010源代码

12、：LU_Decomposition.cpp#include#include#define N 4 /矩阵维数，可自定义static double ANN; /系数矩阵static double BN; /右端项static double YN; /中间项static double XN; /输出static double SN; /选取列主元的比较器int i,j,k; /计数器void main()cout 请输入线性方程组（ai1，ai2，ai3.ain, yi）：endl;for ( i = 0; i N ;i+)for (int j = 0; j Aij;cin Bi;for (k =

13、 0 ; k N; k+)/选列主元int index = k;for (i = k ; i N; i+)double temp = 0;for (int m = 0 ; m k ;m+)temp = temp + Aim*Amk;Si = Aik - temp;if(Sindex Si)index = i;/交换行double temp;for( i = k ; i N ;i+ ) temp = Aindexi; Aindexi = Aki; Aki = temp;temp = Bindex;Bindex = Bk;Bk = temp;/ 构造L、U矩阵for (j = k ; j N; j

14、+)double temp = 0;for (int m = 0 ;m k; m+ )temp = temp + Akm * Amj;Akj = Akj - temp; /先构造U一行的向量for( i = k+1; i N; i+)double temp = 0;for (int m =0 ; m k; m+ )temp = temp + Aim * Amk;Aik = (Aik - temp)/Akk; /再构造L一列的向量/求解LY = BY0 = B0;for (i = 1; i N ; i+)double temp = 0;for (int j =0 ; j = 0 ; i- )do

15、ubletemp = 0;for (int j =i+1 ; j N; j+ )temp = temp + Aij * Xj;Xi = (Yi - temp)/Aii;/打印Xcout 线性方程组的解（X1，X2，X3.Xn）为：endl;for( i = 0; i N ;i+)cout Xi ;运行结果截图：7.简单迭代法（Jacobi迭代）：测试矩阵为：2-100-12-100-12-100-12x1x2x3x4 = 1010 取精度为esp = 0.001，最大迭代次数M = 100源代码：Jacobi.cpp#include#include#define N 4 /矩阵的维数，可按需更

16、改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系数矩阵static double BN=1,0,1,0; /右端项static double YN; /输出比较项static double YN; static double XN; /输出项static double GN; /X = BX + G的G矩阵int i,j,k; /计数器double eps = 0.001;int M = 100;bool distance() /求两输出项的差的范数是否满足精度要求double temp = 0;for (i = 0 ;

17、i eps) return false;elsereturn true; /满足精度要求则结束程序void main() /形成迭代矩阵B，存放到A中for (i = 0 ;i N;i+)if (fabs(Aii) eps)cout 打印失败endl;return;double T = Aii;for ( j = 0 ; j N;j+)Aij = -Aij/T;Aii = 0;Gi = Bi/T;int counter = 0;while (counter M)/迭代for (i = 0;i N; i+)double temp = 0;for (j = 0;jN; j+)temp = temp

18、 + Aij*Yj;Xi = Gi + temp;if (distance()=true)break;else/交换X，Y向量；for( i = 0; i N ;i+)Yi = Xi;counter+;/打印Xcout 迭代次数为：counter次。该线性方程组的解（X1，X2，X3.Xn）为：endl;for( i = 0; i N ;i+)cout Xi ;运行结果截图：8.Seidel迭代：测试矩阵为：2-100-12-100-12-100-12x1x2x3x4 = 1010 取精度为esp = 0.001，最大迭代次数M = 100源代码：Seidel.cpp#include#incl

19、ude#define N 4 /矩阵的维数，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系数矩阵static double BN=1,0,1,0; /右端项static double YN; /输出比较项static double XN; /输出项static double GN; /X = BX + G的G矩阵int i,j,k; /计数器double eps = 0.001;int M = 100;bool distance() /求两输出项的差的范数是否满足精度要求double temp = 0;for

20、(i = 0 ;i eps) return false;elsereturn true; /满足精度要求则结束程序void main() /形成迭代矩阵B，存放到A中for (i = 0 ;i N;i+)if (fabs(Aii) eps)cout 打印失败endl;return;double T = Aii;for ( j = 0 ; j N;j+)Aij = -Aij/T;Aii = 0;Gi = Bi/T;int counter = 0;while (counter M)/迭代for (i = 0;i N; i+)double temp = 0;for (j = 0;jN; j+)tem

21、p = temp + Aij*Xj;Xi = Gi + temp;if (distance()=true)break;else/交换X，Y向量；for( i = 0; i N ;i+)Yi = Xi;counter+;cout 迭代次数为：counter次。该线性方程组的解（X1，X2，X3.Xn）为：endl;for( i = 0; i N ;i+)cout Xi ;9.松弛法（SOR迭代）测试矩阵为： 取松弛系数w = 1.462-100-12-100-12-100-12x1x2x3x4 = 1010 取精度为esp = 0.001，最大迭代次数M = 100源代码：SOR.cpp#inc

22、lude#include#define N 4 /矩阵的维数，可按需更改static double ANN = 2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系数矩阵static double BN = 1,0,1,0; /右端项static double YN; /输出比较项static double XN; /输出项static double GN; /X = BX + G的G矩阵int i,j,k; /计数器double eps = 0.001; /最小容许误差int M = 100; /最大容许比较次数double w = 1.46; /松弛因子bool

23、 distance() /求两输出项的差的范数是否满足精度要求double temp = 0;for (i = 0 ;i eps) return false;elsereturn true; /满足精度要求则结束程序void main() /形成迭代矩阵B，存放到A中for (i = 0 ;i N;i+)if (fabs(Aii) eps)cout 打印失败endl;return;double T = Aii;for ( j = 0 ; j N;j+)Aij = -w * Aij/T;Aii = 1-w;Gi = w * Bi/T;int counter = 0;while (counter M)/迭代for (i = 0;i N; i+)double temp = 0;for (j = 0;jN; j+)temp = temp + Aij*Xj;Xi = Gi + temp;if (distance()=true) break;else/交换X，Y向量；for( i = 0; i N ;i+)Yi = Xi;counter+; cout 迭代次数为：counter次。该线性方程组的解（X1，X2，X3.Xn）为：endl;for( i = 0; i N ;i+)cout Xi ;运行效果截图：24