C++_数值_三次样条插值_自动选取步长梯形法_Romberg求积法_列主元高斯消去法_列主元LU分解法_Jacobi迭

1、数值分析上机实验报告之样条插值1.三次样条插值（初值条件1）：P52.9、给定函数y=fx的函数表和边界条件s''75=0，s''80=0，求三次样条插值函数s(x)，并求f(78.3)的近似值。函数表x757677787980y=fx2.7682.8332.9032.9793.0623.153源代码：yangtiao.cpp#include<iostream.h>#include<math.h>void main()int choice = 0;int n = 2;double xx,*x, *y,*a,*b,*a1,*b1,*h,*m

2、;cout<<"请输入插值节点个数n："<<endl;cin >>n;x = new doublen; y = new doublen;a = 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 ;

3、 i < n ; i+)cin>>xi>>yi;for (int j = 0 ; j < n-1 ; j+)hj = xj+1 - xj;cout<<"请输入待估点xx："<<endl;cin>>xx;cout<<"请选择边界条件："<<endl;cin>>choice;switch(choice)case 1:double temp1,temp2;a0 = 0;an-1 = 1;cout<<"请输入边界条件的两个一阶微商值s

4、'(x1)与s'(xn)："<<endl;cin>>temp1>>temp2;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 ; k<n-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/

5、2;for (int l = 1; l<n ;l+)a1l = -al/(2+(1-al)*a1l-1);b1l = ( bl - (1-al)*b1l-1 )/( 2+ (1-al) * a1l-1 );mn = 0;for (j = n-1 ; j>=0; j-)mj = a1j * mj+1 + b1j;/判别xx所在区间并输出结果cout<<"n插值结果为："for(k = 0 ;k <n-1 ;k+)if (xk <= xx && xk+1 >xx )double output = 0;output = (

6、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;break;delete x;delete y; delete a; delete b; delete a1;delete b1; delete h;

7、 delete m;运行结果截图：2.三次样条插值（初值条件2）：P52.10、给定函数y=fx的函数表和边界条件s'0.25=1，s'0.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+t2dt的近似值。（给定=0.01）源代码：SelfSelLength.cpp#include<iostream.h>#i

8、nclude<math.h>double 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 <= n; k+ )S = S+fun( R_a + (2*k-

9、1)*h/n);R1 = R0/2 + S * h/n;E = fabs(R1 - R0);n = 2*n;while(E > 3*e);return R1;void main()double a,b,e;cout << "请依次输入待求积分函数的下界a、上界b 及 精度要求e:" << endl;cin >> a >> b >> e;cout <<"自动选取步长梯形法可求得积分:"<< SelfSelLength (a,b,e)<<endl;运行结果截

10、图：4.Romberg求积法：P978、使用Romberg求积法计算积分19xdx的近似值。（给定=0.01，且取n=1）源代码：Romberg.cpp#include<iostream.h>#include<math.h>static double Tri128128;double fun(double a)return sqrt(a);double Romberg(double 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 /每当误差值不符合要

11、求时，计算下一行的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 <= k ; m +)Trimk-m = ( pow(4,m)*Trim-1k-m+1 - Trim-1k-m ) / (pow(4,m) - 1);E = fabs(Trik0 -

12、 Trik-10);while(E > e);return Trik0;void main()double a,b,e;cout << "请依次输入待求积分函数的下界a、上界b 及 精度要求e:" << endl;cin >> a >> b >> e;cout <<"按Romberg求积法可求得积分:"<< Romberg (a,b,e)<<endl;运行结果截图：5.列主元高斯消去法：测试矩阵为：2-100-12-100-12-100-12×

13、x1x2x3x4 = 1010源代码：Guess_Elimination.cpp/列主元高斯消去法#include<iostream.h>#include<math.h>#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 XN; int i,j,k; /计数器void main()for(k = 0; k < N-1 ;k+)/选取最大主元int index

14、 = 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 = temp;for(i = k+1; i<N; i+)double T = Aik/Akk;Bi = Bi - T * Bk;for ( j = k+1 ; j < N ; j+ )Aij = Aij - T * Ak

15、j;XN-1 = BN-1/AN-1N-1;for (i = N-2; i >=0 ; i-)double Temp = 0;for (int j = i+1; j<N ;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

16、-12-100-12×x1x2x3x4 = 1010源代码：LU_Decomposition.cpp#include<iostream.h>#include<math.h>#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

17、，ai2，ai3.ain, yi）："<<endl;for ( i = 0; i < N ;i+)for (int j = 0; j< N ;j+ )cin >> Aij;cin >>Bi;for (k = 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)ind

18、ex = 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+)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

19、 = 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 < i; j+ )temp = temp + Aij * Yj;Yi = Bi - temp;/求解UX = YXN-1 = YN-1/AN-1N-1;for (i = N-2 ; i >= 0 ; i- )doubletemp = 0;for (

20、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-12×x1x2x3x4 = 1010 取精度为esp = 0.001，最大迭代次数M =

21、100源代码：Jacobi.cpp#include<iostream.h>#include<math.h>#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 YN; static double XN; /输出项static double GN; /X = BX' + G的G矩阵int i,j,k; /计数器d

22、ouble eps = 0.001;int M = 100;bool distance() /求两输出项的差的范数是否满足精度要求double temp = 0;for (i = 0 ;i< N ; i+)temp = temp + fabs(Xi - Yi);if (temp > eps) return false;elsereturn true; /满足精度要求则结束程序void main() /形成迭代矩阵B，存放到A中for (i = 0 ;i< N;i+)if (fabs(Aii) < eps)cout <<"打印失败"<

23、<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;j<N; j+)temp = temp + Aij*Yj;Xi = Gi + temp;if (distance()=true)break;else/交换X，Y向量；for( i = 0; i < N ;i+)Yi = Xi;co

24、unter+;/打印Xcout << "迭代次数为："<<counter<<"次。该线性方程组的解（X1，X2，X3.Xn）为："<<endl;for( i = 0; i < N ;i+)cout << Xi <<" "运行结果截图：8.Seidel迭代：测试矩阵为：2-100-12-100-12-100-12×x1x2x3x4 = 1010 取精度为esp = 0.001，最大迭代次数M = 100源代码：Seidel.cpp#include&l

25、t;iostream.h>#include<math.h>#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()

26、 /求两输出项的差的范数是否满足精度要求double temp = 0;for (i = 0 ;i< N ; i+)temp = temp + fabs(Xi - Yi);if (temp > 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 ;

27、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;j<N; 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 << "迭代次数为："<&l

28、t;counter<<"次。该线性方程组的解（X1，X2，X3.Xn）为："<<endl;for( i = 0; i < N ;i+)cout << Xi <<" "9.松弛法（SOR迭代）测试矩阵为： 取松弛系数w = 1.462-100-12-100-12-100-12×x1x2x3x4 = 1010 取精度为esp = 0.001，最大迭代次数M = 100源代码：SOR.cpp#include<iostream.h>#include<math.h>#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