实验一实验项目：共轭梯度法求解对称正定的线性方程组 实验容：用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k xx --∞-<时停止计算。

编制程序：储存m 文件function [x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; elsebeta=rho/rho1; p=r+beta*p; end w=A*p;alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho; rho=r'*r; end运行程序：clear,clcA=[2 -1 0 0;-1 3 -1 0;0 -1 4 -1;0 0 -1 5]; b=[3 -2 1 5]'; [x,k]=CGmethod(A,b)运行结果： x =1.1176 0.2353 0.5882 1.1177(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵，b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。

编制程序：储存m 文件function [x,k]=CGmethod_1(A,b) n=length(A);x(1:n,1)=0;r=b-A*x;r1=r; k=0;while norm(r1,1)>10^(-7)&k<10^4 k=k+1; if k==1 p=r; elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end运行程序：clear,clc n=1000; A=hilb(n); b=sum(A')'; [x,k]=CGmethod(A,b)实验二1、 实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。

实验容：计算下列定积分(1) dx x x x ⎰⎪⎪⎭⎫ ⎝⎛+-202610 (2)⎰10dx x x (3) ⎰20051dx x 实验要求：（1）分别用复化Simpson 公式、自适应复化梯形公式计算要求绝对误差限为71021-⨯=ε，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；（2）分析比较计算结果。

程序：syms xf=x^6/10-x^2+x %定义函数f （x ） n=input('输入所求导数阶数:')f2=diff(f,x,n) %求f(x)的n 阶导数（1）复化梯形 clcsyms x%定义自变量xf=inline('x^6/10-x^2+x','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('3*x^4 - 2','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。

f3='-(3*x^4 - 2)'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=0.5*10^(-7) %精度要求值a=0 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求连加和xk=a+k*hTn1=Tn1+f(xk)endTn=h/2*((f(a)+2*Tn1+f(b)))z=exp(2)R=Tn-z %求已知值与计算值的差stem(xk,Tn1);fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn)fprintf('等分数 n=')disp(n) %输出等分数fprintf('已知值与计算值的误差 R=')disp(R)（2）复化Simpsonclcclearsyms x%定义自变量xf=inline('x^6/10-x^2+x','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('36*x^2','x') %定义f(x)的四阶导数，输入程序1里求出的f2即可f3='-(36*x^2)'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10^(-8) %精度要求值a=0 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endh=(b-a)/n %求h Sn1=0 Sn2=0for k=0:n-1 %求两组连加和 xk=a+k*h xk1=xk+h/2 Sn1=Sn1+f(xk1) Sn2=Sn2+f(xk) endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b)) %因Sn2多加了k=0时的值，故减去f(a) z=exp(2)R=Sn-z %求已知值与计算值的差 fprintf('用Simpson 公式计算的结果 Sn=') disp(Sn)fprintf('等分数 n=')disp(n) fprintf('已知值与计算值的误差 R=') disp(R)结果：dx x x x ⎰⎪⎪⎭⎫ ⎝⎛+-202610 用复化梯形算法计算的结果 Tn= 1.6674等分数 n= 24764已知值与计算值的误差 R= -6.3976用Simpson 公式计算的结果 Sn= 1.0434等分数 n= 76已知值与计算值的误差 R= -6.0216⎰1dx x x用复化梯形算法计算的结果 Tn= 0.8985等分数 n= 1119已知值与计算值的误差 R= -6.1665用Simpson 公式计算的结果 Sn= 0.9406已知值与计算值的误差 R= -6.1245⎰20051dxx用复化梯形算法计算的结果 Tn= 23.2353 等分数 n= 1000000已知值与计算值的误差 R= 16.1704 用Simpson 公式计算的结果 Sn= 23.2617 等分数 n= 10647已知值与计算值的误差 R= 16.1969分析：在处理问题时，复化Simpson 要比复化梯度计算速度要快很多。

2、实验目的：高斯数值积分方法用于积分方程求解。

实验容：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。

而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。

在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。

对第二类Fredholm 积分方程b t a t f ds s y s t k t y ba≤≤+=⎰),()(),()(首先将积分区间[a,b]等分成n 份，在每个子区间上离散方程中的积分就得到线性代数方程组。

实验要求：分别使用如下方法，离散积分方程中的积分1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。

求解如下的积分方程t te ds s y e e t y --=⎰10)(12)(，方程的准确解为t e ， 并比较各算法的优劣。

程序结果：当迭代次数复化梯形方法的平均误差err=0.247复化辛甫森方法的平均误差err=0.00116 复化高斯方法的平均误差err=0.0008复化梯形方法的平均误差err=0.00116复化辛甫森方法和复化高斯方法的平均误差err=0可以看出，复化高斯方法得到的结果精度最高，复化辛普森方法比复化梯形方法的精度要高，程序：clear,clca=0;b=1;n=5;figurefun1(a,b,n);hold ona=0;b=1;n=5;figurefun2(a,b,n);hold onfigurefun3(a,b,n);编制m函数：function y=Kfun(t,s)y=2/(exp(1)-1)*exp(t);function y=ffun(t)y=-exp(t);function y=Fexc(t)%精确解y=exp(t);function [x,w]=fhgauss(a,b,n)h=(b-a)/n;x1=a:h:b;x=zeros(1,2*n);w=x;for i=1:n[x(2*i-1:2*i),w(2*i-1:2*i)]=GaussLegendre(x1(i),x1(i+1),2);endx=a:h/2:b;w=x;w(1)=h/6;w(2*n+1)=h/6;for i=1:nw(2*i)=2/3*h;if i<nw(2*i+1)=h/3;endendfunction [x,w]=fhtrapz(a,b,n)h=(b-a)/n;x=a:h:b;for i=1:n+1if i==1||i==n+1w(i)=h/2;elsew(i)=h;endendfunction [x,w]=GaussLegendre(a,b,n) i=1:n-1;c=i./sqrt(4*i.^2-1);CM=diag(c,1) + diag(c,-1);[V L]=eig(CM);[x ind]=sort(diag(L));V=V(:,ind)';w=2*V(:,1).^2;x=x*(b-a)/2+(b+a)/2;w=w*(b-a)/2;function fun1(a,b,n)[x1,w1]=fhtrapz(a,b,n);n1=4;n=n+1;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i));for j=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j)); endendA=eye(n)-A;y1=(A\B)';yN=x;for i=1:n1+1yN(i)=ffun(x(i));for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j); endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');function fun2(a,b,n)[x1,w1]=fhsimpson(a,b,n);n=2*n+1;n1=50;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i));endfor i=1:nfor j=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j)); endendA=eye(n)-A;y1=(A\B)';yN=x;for i=1:n1+1yN(i)=ffun(x(i));for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j);endfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');function fun3(a,b,n)[x1,w1]=fhgauss(a,b,n);n=2*n;n1=4;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i));endfor i=1:nfor j=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j)); endendA=eye(n)-A;y1=(A\B)';yN=x;for i=1:n1+1yN(i)=ffun(x(i));for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j); endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');图一图二图三实验三1、对常微分方程初值问题2cos (0)1y y x y '=-+⎧⎨=⎩ (0)x π<< 分别使用Euler 显示方法、改进的Euler 方法和经典RK 法和四阶Adams 预测-校正算法，求解常微分方程数值解，并与其精确解cos sin y x x =+进行作图比较。