浙大计算方法上机报告

1、学号：3100300038 姓名:作业1:用列主元高斯消去法和列主元三角分解法解P227页第3题1.列主元高斯消去法目的：用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现akk =0情况，这时必须进行带行交换的高斯消去法。但在实际计算中即使akk0但其绝对值很小时，用akk作除数，会导致中间结果矩阵A(k)元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠。 列主元高斯消去法可以难过一般高斯法的这些缺点。一、列主元高斯消去法解方程的Matlab程序如下：function a=columneli(a)%对矩阵a进行列主元消去m n=size(a);%求取 a 的行数 m禾D歹U

2、数 nfor i=1:m-1 ImaxEle ,p os=max(abs(a(i:e nd,i);maxRow=po s+i-1;if a(maxRow,i)=0disp(矩阵为奇异矩阵)returnendif(maxRow=i)temp=a(maxRow,:);a(maxRow,:)=a(i,:);a(i,:)=te mp; end for j=i+1:ma(j,i)=a(j,i)/a(i,i);for k=i+1: n在每次变换前寻找绝对值最大主所在列maxRow对于非奇异矩阵，在程序中给予提示，结束程序与列主元绝对值最大的行进行行交换求取第j列主元对第j列主元进行行变换利用列主元消去的结

3、果求方程的解,a为方程组的增广矩阵a(j,k)=a(j,k)-a(j,i)*a(i,k); %endendendfunction x=elisolve(a) % a=colu mn eli(a);m, n=size(a);x=zeros(m,1);for i=m:-1:1x(i)=(a(i, n)-a(i,i:m)*x(i:m)/a(i,i); end列主元高斯消去法解P227页第3题:01-1L2答案:仁24x112jLx3j L3jx2x=7/6 -1/3 1/2程序流程图如下:酉用圉ft时亘诅订化 寸-恵冷去 fi-col umitinl i (rt)条件图框里面没有条件式，因为i是从m

4、到1,所以i=1后运行下面的命令，直接输出结杲。M 亘划 uWinpIN*)r*ilT,flKRQT*rr J f a :) -a (i j : 11 3Li.L, ： tcJiipv u交議lT，iixrw*rF syni(ans)arts：/6-1/31/2运行结果：x=7/6 -1/3 1/2 上机结果和正确答案完全相同。2.列主元三角分解法一，列主元三角分解法解方程的Matlab程序如下：function l,u,p=tridecom(a)%对矩阵 a 进行三角分解m n=size(a); % 求取a的行数 m和列数n if(m=n)dis p(矩阵的行数与列数不相等，无法进行三角分解

5、returnendl=zeros (n );u=zeros (n); p=eye( n);%p *a=l*umaxEle,goodRow=max(abs(a(:,1);行交换temp=a(1,:);a(1,:)=a(goodRow,:);a(goodRow,:)=te mp; temp=p (1,:); p(1,:)=p( goodRow,:); p( goodRow,:)=te mp; % u(1,:)=a(1,:);l(1:e nd,1)=a(1:e nd,1)/u(1,1);for r=2: ntempU rr=zeros (n ,1); for k=r: ntemp Urr(k)=a(

6、k,r)-l(k,1:r-1)*u(1:r-1,r);endmaxEle,goodRow=max(abs(te mpU rr);temp=a(r,:);a(r,:)=a(goodRow,:);a(goodRow,:)=te mp; temp=l(r,:);l(r,:)=l(goodRow,:);l(goodRow,:)=tem p; temp=p( r,:); p( r,:)=p( goodRow,:); p( goodRow,:)=te mp; for i=r: nu(r,i)=a(r,i)-l(r,1:r-1)*u(1:r-1,i);endfor i=r: nl(i,r)=(a(i,r)-

7、l(i,1:r-1)*u(1:M,r)/u(r,r);endendfunction x=lusolve(a) %m, n=size(a);l,u,p=tndecom(a(:,1:m); % x=p *a(:, n);for i=2:m用三角分解法解线性方程组Ax=b , a为方程组的增广矩阵A|b对A进行三角分解 A=p*l*u对a的最后一列进行与方程组系数三角分解时相同的行变换求出(|A-1)*p*b求出方程组的解 x=(UA-1)*(LA-1)* p*bx(i)=x(i)-l(i,1:i-1)*x(1:i-1); %endx(m)=x(m)/u(m,m);for i=m-1:-1:1x(i

8、)=(x(i)-u(i,i+1:m)*x(i+1:m)/u(i,i);%end、程序流程图如下:t= 1V f III. Ill 、件堆广虬PT .|rH治|ir, -size m；L| Ur p -t ride CO- I a , 13； t ill ；;对肌出冇的解為-p -y 一x-. a FL!;林二的轻F；列i宜吁程阻癢数：肃汀册比ttl诃弗n建换X III i )1= X U 11 C i 11: i 三 1】* (1: i 三 1】；ft ：iN 1 -11 *p*cX t r ； X ijii j j L U； f ir.1 ixH) - IIIX !1 i - U | 1,

9、i f L : 111) x il r i ；m) J A; !1, 1)r.nalYESMO1 J Hr pTft I !埠忧石-【Al bT 4绘 T输出*卯阵的行s与刊Brf*租導00p E严nlIhwhELr -；-:-乐gsi (4 I: J-lrO:indxE 1 f y11 t.j Jt i -tai pl i i ml, t J a El, 1 Jll :严厂叫.J = ?r：1, 1. ) . E： ( 12 : nk r : Ttcp-j-rr Elk )L(k, 121-1! *,1 U r 13 (包 t i . )-亠(i ,二 3 上 t I . I. : r-1

10、r ) .1 /i； r, T ；,;(丿；K qpnHS4r-Wj1 j r Ij UL p =t ri且亡匸4血 D 31.1 t i : 21 Z)1.0000J.OOOQ00- 5000-0- 5000 I- tJOfliO可见，p*a=l*u作业2；编写Matlab程序，用最小二乘法拟合多项式曲线（解P89页第1题）目的和意义；根据观测或实验得到一系列的数据，确定了与自变量的某些点相应的函数值。当函数比较复杂或根本无法写出解析式时，往往根据观测数据构造一个 适当的简单的函数近似地代替要寻求的函数。利用最小二乘法拟合多项式来描述其函数。一、用最小二乘法拟合多项式曲线的Matlab程序如

11、下:fun cti on a=po lyercheng(x,y,n)x=x;y=y;m=size(x,1);c=on es(m,1+ n);for i=1: nc(:,i+1)=x.Ai;end.A=c*c;b=c*y;a=(Ab); Ia=fl ip lr(a);|c=pol yval(a,x)-y;d=max(abs(c);%fn,d);fprintfC最大偏差：d=sqrt(c*c);fprintfC均方误差：b=li nsp ace( min( x),max(x),500); c=pol yval(a,b);plot(x,y,g*,b,c,r);end、程序的流程图如下：_ a=po

12、lyercheng(x,y,n)i开始输入一组数据(x,y)和多项式次数nx=x;y=y;m=size(x,1);生成中间矩阵Cc=on es(m, n+1);I4一 i=1：nC(:,i+1)=x.Ai;IA=c*c;b=c*y;Ja=(Ab);a=fli plr(a);生成偏差阵列解方程组得多项式系数向量fpr in tf(d=max(abs(c);最大偏差：%fn,d);1fprintf(d=sqrt(c*c);均方误差：%f,d);1输出a结束程序c=pol yval(a,x)-y;三、最小二乘法拟合多项式曲线解 P89页第1题: 已知实验数据如下：x2468Y2112840用最小二乘

13、法求一次和二次拟合多项式，分别算出均方误差和最大偏差。 结果：二次拟合多项式最大偏差：1.950000均方误差：2.906888拟合多项式曲线：y=0.1875x2+4.6750x-8.7500一次拟合多项式最大偏差：2.700000均方误差：3.271085拟合多项式曲线：y=6.5500x-12.5000四、解P89页第1题在Matlab中的截图如下:2. 7000003.2；10S5 a=polyer cheng y, 1) 最大偏羞； 均方误差；6,5500-12.5000y=6.5500x-12.50001.9500002, 906S38 a=polyercheng (x, y, 2

14、、 最大偏差： 均方误差：Q18754.6750乩 7500y=0.1875x2+4.6750x-8.7500其中a为拟合多项式的系数阵列。作业3 :用龙贝格算法计算，fdx使截断误差不超过0.5e-6。0 x目的和意义；在一元函数的积分学中，我们已经熟知，若函数f(x)通过最小二乘法拟合多项式得到的多项式接近与实验值(代入后)。续且其原函数为F(x),则可用牛顿一莱布尼兹公在区间a, b 上连 bf(x)dx=F(a)-F(b)来求定积分。牛顿一莱布尼兹公式虽然在理论上或在解决实际问题中都起 了很大的作用，但它并不能完全解决定积分的计算问题。当被积函数f(x)在区间a,b 上连续时,要使得复

15、合梯形公式比较精确 地代替定积分bJ f (x)dx = F(a) - F(b)可将分点(即基点)加密,也就是将区间a,b a细分,然后利用复合梯形公式求积将收敛缓慢的梯形值序列加工成收敛迅速的Sn,Cn,Rn，这种加速方法称为龙贝格算法。它可以把复杂的积分比较精确的求出来。、龙贝格算法Matlab程序如下：function z=romberg(ff,a,b,e psi,k) %if nargin=4k=16;end%t=zeros(1,k+1);format Io ng;fun ction y=f(x)%嵌套函数，用来判断积分区域内的奇点y=ff(x); %若 ff 参数为 (x)sin(x

16、)/x,则有 y=sin(x)/xif isnan(y)=true|is inf(y )=truedis p(请确保所给函数在积分区域内无奇点);endendt(1)=(b-a)/2*(f(a)+f(b);for l=1:ksum=0;for i=1:(2A(l-1)sum=sum+f(a+(2*i-1)*(b-a)/(2A|);end用龙贝格算法计算积分如果只有4个传入参数，则将二分次数 k设置为默认值16t(l+1)=t(l)/2+(b-a)/(2A|)*sum;endfor m=1:ktemp=t;for i=m+1:k+1%将新序列放入 t(m+1)-t(k+1) 中t(i)=(4m)*te mp(i)-te mp (卜1)/(4人口-1);|end计算梯形序列if abs(t(m+1)-t(m) rojnbergsin(xl /x, epSj 1, 0. 5e-6)截断误差为S. S32355e-008ans =0-9