数值分析实验报告.docx

1、数值分析实验报告 数值分析实验报告系 别电子信息系专 业计算机科学与技术班级学号姓 名指导教师2011年 6 月 20日实验一实验题目：编写一个拉格朗日插值函数，对不多于9个点的插值节点都可以求出插值函数，任意给定输入x值都可以求出y值。程序代码：#include #include #include using namespace std;#define N 100 void lagrange() int n,k,m,q=1; float xN,yN,xx,yyy1,yyy2,yy1,yy2,yy3; coutn; for(k=0;k=n-1;k+) cout请输入Xkxk; cout请输入Y

2、kyk; system(cls); cout则Xi与Yi表格如下：endl; coutXi ;for(k=0;k=n-1;k+)coutsetiosflags(ios:left)setw(10)xk; coutendl; coutYi ;for(k=0;k=n-1;k+)coutsetiosflags(ios:left)setw(10)yk; coutendl; while(q) coutxx; while(xxxk-1|xxx0) coutxx; for(k=0;k=n-1;k+) if(xxxk) m=k-1; k=n-1; yyy1=ym*(xx-xm+1)/(xm-xm+1)+ym+1

3、*(xx-xm)/(xm+1-xm); cout则拉格朗日分段线性插值为：yyy1endl; for(k=0;k=n-1;k+) if(xx(xm+1-xx)m=m+1; else m=m; yy1=ym-1*(xx-xm)*(xx-xm+1)/(xm-1-xm)*(xm-1-xm+1); yy2=ym*(xx-xm-1)*(xx-xm+1)/(xm-xm-1)*(xm-xm+1); yy3=ym+1*(xx-xm-1)*(xx-xm)/(xm+1-xm-1)*(xm+1-xm); yyy2=yy1+yy2+yy3; cout则拉格朗日分段二次插值为：yyy2endl; coutq; syst

4、em(cls);void main() lagrange();实验结果：实验二实验题目：用不同的方法计算积分 取不同的步长h，分别用复合梯形公式及复合辛普森公式计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善。程序代码：复合梯形公式：function I = T_quad(x,y)n=length(x); m=length(y);if n=merror(The lengths of X and Y must be equal);return;endh=(x(n)-x(1)/(n-1);a=1 2*ones(1,n-2) 1;I=h/

5、2 * sum(a.*y);复合辛普森公式：function I = S_quad(x,y)n = length(x);m = length(y);if n=merror(The lengths of X and Y must be equal);return;endif rem(n-1,2)=0I = T_quad(x,y);return;endN = (n-1)/2;h=(x(n)-x(1)/N;a = zeros(1,n);for k =1:Na(2*k - 1) = a(2*k-1)+1;a(2*k) = a(2*k)+4;a(2*k + 1) = a(2*k+1)+1;endI =

6、h/6 * sum(a.*y);测试数据：复合梯形测试数据： x=0.00001:0.0001:0.99999 y=sqrt(x).*log (x)I=T_quad(x,y)复合辛普森测试数据： x=0.00001:0.0001:0.99999 y=sqrt(x).*log (x)I=S_quad(x,y)实验结果：复合梯形实验结果： x=0.00001:0.0001:0.99999x = Columns 1 through 8 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 Columns 9 through 16 Column

7、s 9993 through 10000 0.9992 0.9993 0.9994 0.9995 0.9996 0.9997 0.9998 0.9999 y=sqrt(x).*log (x)y = Columns 1 through 8 -0.0364 -0.0956 -0.1227 -0.1422 -0.1579 -0.1712 -0.1828 -0.1932 Columns 9993 through 10000 -0.0008 -0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 I=T_quad(x,y)I = -0.4444I

8、= -0.4444复合辛普森实验结果： x=0.00001:0.0001:0.99999x = Columns 1 through 8 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007Columns 9993 through 10000 0.9992 0.9993 0.9994 0.9995 0.9996 0.9997 0.9998 0.9999 y=sqrt(x).*log (x)y = Columns 1 through 8 -0.0364 -0.0956 -0.1227 -0.1422 -0.1579 -0.1712 -0.1

9、828 -0.1932 Columns 9993 through 10000 -0.0008 -0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 I=S_quad(x,y)I = -0.4444I = -0.4444实验三实验题目：用LU分解和列主元消去法解线性方程组输出Ax=b中系数A=LU分解的矩阵L和U，解向量x及detA；列主元法的行交换次序，解向量x及detA；比较两种方法所得的结果。程序代码：/解线性方程组#include#include#include/-全局变量定义区const int Number=15; /方程最大

10、个数double aNumberNumber,bNumber,copy_aNumberNumber,copy_bNumber; /系数行列式int A_yNumber; /a中随着横坐标增加列坐标的排列顺序,如a00,a12,a21.则A_y=0,2,1.;int lenth,copy_lenth; /方程的个数double a_sum; /计算行列式的值char * x; /未知量a,b,c的载体/-函数声明区void input(); /输入方程组void print_menu(); /打印主菜单int choose (); /输入选择void cramer(); /Cramer算法解方程

11、组void gauss_row(); /Gauss列主元解方程组void Doolittle(); /用Doolittle算法解方程组int Doolittle_check(double aNumber,double bNumber); /判断是否行列式0,若是,调整为顺序主子式全0void xiaoqu_u_l(); /将行列式Doolittle分解void calculate_u_l(); /计算Doolittle结果 double & calculate_A(int n,int m); /计算行列式double quanpailie_A(); /根据列坐标的排列计算的值,如A_y=0,2

12、,1,得sum=a0 A_y0 * a1 A_y1 * a2 A_y2 =a00*a12*a21;void exchange(int m,int i); /交换A_ym,A_yivoid exchange_lie(int j); /交换aj与b;void exchange_hang(int m,int n); /分别交换a和b中的m与n两行void gauss_row_xiaoqu(); /Gauss列主元消去法void gauss_calculate(); /根据Gauss消去法结果计算未知量的值void exchange_a_lie(int m,int n); /交换a中的m和n列void

13、 exchange_x(int m,int n); /交换x中的xm和xnvoid recovery(); /恢复数据/主函数void main() int flag=1; input(); /输入方程 while(flag) print_menu(); /打印主菜单 flag=choose(); /选择解答方式 /函数定义区void print_menu() system(cls); cout-方程系数和常数矩阵表示如下:n; for(int j=0;jlenth;j+) cout系数j+1 ; coutt常数; coutendl; for(int i=0;ilenth;i+) for(j=

14、0;jlenth;j+) coutsetw(8)setiosflags(ios:left)aij; couttbiendl; cout-请选择方程解答的方案-; coutn 1. Gauss列主元消去法; coutn 2. Doolittle分解法;coutn 3. 退出; coutn 输入你的选择:;void input() int i,j; coutlenth; if(lenthNumber) coutIt is too big.n; return; x=new charlenth; for(i=0;ilenth;i+) xi=a+i;/输入方程矩阵 /提示如何输入 cout=n; cou

15、t请在每个方程里输入lenth系数和一个常数:n; cout例：n方程:a; for(i=1;ilenth;i+) cout+i+1xi; cout=10n; cout应输入:; for(i=0;ilenth;i+) couti+1 ; cout10n; cout=n; /输入每个方程 for(i=0;ilenth;i+) cout输入方程i+1:; for(j=0;jaij; cinbi; /备份数据 for(i=0;ilenth;i+) for(j=0;jlenth;j+) copy_aij=aij; for(i=0;ichoice; switch(choice) case 1:crame

16、r();break; case 2:gauss_row();break; case 3:guass_all();break; case 4:Doolittle();break; case 5:return 0; default:cout输入错误,请重新输入:; choose(); break; coutch; if(ch=n|ch=N) return 0; recovery(); coutnnn; return 1;/高斯列主元排列求解方程void gauss_row() int i,j; gauss_row_xiaoqu(); /用高斯列主元消区法将系数矩阵变成一个上三角矩阵 for(i=0

17、;ilenth;i+) for(j=0;jlenth;j+) coutsetw(10)setprecision(5)aij; coutsetw(10)biendl; if(alenth-1lenth-1!=0) cout系数行列式不为零,方程有唯一的解：n; gauss_calculate(); for(i=0;ilenth;i+) /输出结果 coutxi=bin; else cout系数行列式等于零,方程没有唯一的解.n;void gauss_row_xiaoqu() /高斯列主元消去法 int i,j,k,maxi;double lik; cout用Gauss列主元消去法结果如下:n;

18、for(k=0;klenth-1;k+) j=k; for(maxi=i=k;iamaxij) maxi=i; if(maxi!=k) exchange_hang(k,maxi);/ for(i=k+1;ilenth;i+) lik=aik/akk; for(j=k;jlenth;j+) aij=aij-akj*lik; bi=bi-bk*lik; void Doolittle() /Doolittle消去法计算方程组 double temp_aNumberNumber,temp_bNumber;int i,j,flag; for(i=0;ilenth;i+) for(j=0;jlenth;j

19、+) temp_aij=aij; flag=Doolittle_check(temp_a,temp_b); if(flag=0) coutn行列式为零.无法用Doolittle求解.; xiaoqu_u_l(); calculate_u_l(); cout用Doolittle方法求得结果如下:n; for(i=0;ilenth;i+) /输出结果 for(j=0;xj!=a+i&jlenth;j+); coutxj=bjendl; void calculate_u_l() /计算Doolittle结果 int i,j;double sum_ax=0; for(i=0;ilenth;i+) fo

20、r(j=0,sum_ax=0;j=0;i-) for(j=i+1,sum_ax=0;jlenth;j+) sum_ax+=aij*bj; bi=(bi-sum_ax)/aii; void xiaoqu_u_l() /将行列式按Doolittle分解 int i,j,n,k;double temp; for(i=1,j=0;ilenth;i+) aij=aij/a00; for(n=1;nlenth;n+) /求第n+1层的上三角矩阵部分即U for(j=n;jlenth;j+) for(k=0,temp=0;kn;k+) temp+=ank*akj; anj-=temp; for(i=n+1;

21、ilenth;i+) /求第n+1层的下三角矩阵部分即L for(k=0,temp=0;kn;k+) temp+=aik*akn; ain=(ain-temp)/ann; int Doolittle_check(double temp_aNumber,double temp_bNumber) /若行列式不为零,将系数矩阵调整为顺序主子式大于零 int i,j,k,maxi;double lik,temp; for(k=0;klenth-1;k+) j=k; for(maxi=i=k;itemp_amaxij) maxi=i; if(maxi!=k) exchange_hang(k,maxi);

22、 for(j=0;jlenth;j+) temp=temp_akj; temp_akj=temp_amaxij; temp_amaxij=temp; for(i=k+1;ilenth;i+) lik=temp_aik/temp_akk; for(j=k;jlenth;j+) temp_aij=temp_aij-temp_akj*lik; temp_bi=temp_bi-temp_bk*lik; if(temp_alenth-1lenth-1=0) return 0; return 1;void exchange_hang(int m,int n) /交换a中和b两行 int j; double

23、 temp; for(j=0;jlenth;j+) temp=amj; amj=anj; anj=temp; temp=bm; bm=bn; bn=temp;void exchange(int m,int i) /交换A_ym,A_yi int temp; temp=A_ym; A_ym=A_yi; A_yi=temp;void exchange_lie(int j) /交换未知量b和第i列 double temp;int i; for(i=0;ilenth;i+) temp=aij; aij=bi; bi=temp; void exchange_a_lie(int m,int n) /交换a中的两列 double temp;int i; for(i=0;ilenth;i+) temp=aim; aim=ain; ain=temp; void exchange_x(int m,int n) /交换未知量xm与xn char temp; temp=xm; xm=xn; xn=temp;void recovery() /用其中一种方法求解后恢复数据以便用其他方法求解 for(int i=0;ilenth;i+) for(int j=0;jlenth;j+) ai