当前位置：文档之家 › 浙大计算方法上机报告.

浙大计算方法上机报告.

  • 格式：doc
  • 大小：1.67 MB
  • 文档页数：14

学号：3100300038 姓名： 专业：作业1：用列主元高斯消去法和列主元三角分解法解P227页第3题1. 列主元高斯消去法目的：用高斯消去法解Ax=b 时，其中设A 为非奇异矩阵，可能出现a kk =0 情况，这时必须进行带行交换的高斯消去法。

但在实际计算中即使0)( k kk a 但其绝对值很小时，用)(k kk a 作除数，会导致中间结果矩阵A(k)元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠。

列主元高斯消去法可以难过一般高斯法的这些缺点。

一、列主元高斯消去法解方程的Matlab 程序如下：function a=columneli(a) %对矩阵a 进行列主元消去 [m n]=size(a); %求取a 的行数m 和列数n for i=1:m-1[maxEle,pos]=max(abs(a(i:end,i)));maxRow=pos+i-1; %在每次变换前寻找绝对值最大主所在列maxRow if a(maxRow,i)==0disp('矩阵为奇异矩阵')return %对于非奇异矩阵，在程序中给予提示,结束程序 end if(maxRow~=i)temp=a(maxRow,:);a(maxRow,:)=a(i,:);a(i,:)=temp;end %与列主元绝对值最大的行进行行交换 for j=i+1:ma(j,i)=a(j,i)/a(i,i); %求取第j 列主元 for k=i+1:na(j,k)=a(j,k)-a(j,i)*a(i,k); %对第j 列主元进行行变换 end end endfunction x=elisolve(a) %利用列主元消去的结果求方程的解,a 为方程组的增广矩阵 a=columneli(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题：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-321321212111430x x x 答案：x=[7/6 -1/3 1/2]T三、程序流程图如下：条件图框里面没有条件式，因为i 是从m 到1，所以i=1后运行下面的命令，直接输出结果。

这里也是一样的，直接运行下面的命令。

四、Matlab中的命令截图如下：运行结果：x=[7/6 -1/3 1/2]T上机结果和正确答案完全相同。

2. 列主元三角分解法一，列主元三角分解法解方程的Matlab程序如下：function [l,u,p]=tridecom(a) %对矩阵a进行三角分解[m n]=size(a); %求取a的行数m和列数nif(m~=n)disp('矩阵的行数与列数不相等，无法进行三角分解')returnendl=zeros(n);u=zeros(n);p=eye(n); %p*a=l*u[maxEle,goodRow]=max(abs(a(:,1)));temp=a(1,:);a(1,:)=a(goodRow,:);a(goodRow,:)=temp;temp=p(1,:);p(1,:)=p(goodRow,:);p(goodRow,:)=temp; %行交换u(1,:)=a(1,:);l(1:end,1)=a(1:end,1)/u(1,1);for r=2:ntempUrr=zeros(n,1);for k=r:ntempUrr(k)=a(k,r)-l(k,1:r-1)*u(1:r-1,r);end[maxEle,goodRow]=max(abs(tempUrr));temp=a(r,:);a(r,:)=a(goodRow,:);a(goodRow,:)=temp;temp=l(r,:);l(r,:)=l(goodRow,:);l(goodRow,:)=temp;temp=p(r,:);p(r,:)=p(goodRow,:);p(goodRow,:)=temp;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)-l(i,1:r-1)*u(1:r-1,r))/u(r,r);endendfunction x=lusolve(a) %用三角分解法解线性方程组Ax=b，a为方程组的增广矩阵[A|b] [m,n]=size(a);[l,u,p]=tridecom(a(:,1:m)); %对A进行三角分解A=p*l*ux=p*a(:,n); %对a的最后一列进行与方程组系数三角分解时相同的行变换for i=2:mx(i)=x(i)-l(i,1:i-1)*x(1:i-1); %求出(l^-1)*p*bendx(m)=x(m)/u(m,m);for i=m-1:-1:1x(i)=(x(i)-u(i,i+1:m)*x(i+1:m))/u(i,i);%求出方程组的解x=(U^-1)*(L^-1)*p*b end二、程序流程图如下：二、解P227页第3题在Matlab中的截图如下：可见，p*a=l*u作业2；编写Matlab程序，用最小二乘法拟合多项式曲线（解P89页第1题）目的和意义；根据观测或实验得到一系列的数据，确定了与自变量的某些点相应的函数值。

当函数比较复杂或根本无法写出解析式时，往往根据观测数据构造一个适当的简单的函数近似地代替要寻求的函数。

利用最小二乘法拟合多项式来描述其函数。

一、用最小二乘法拟合多项式曲线的Matlab程序如下：function a=polyercheng(x,y,n)x=x';y=y';m=size(x,1);c=ones(m,1+n);for i=1:nc(:,i+1)=x.^i;endA=c'*c;b=c'*y;a=(A\b)';a=fliplr(a);c=polyval(a,x)-y;d=max(abs(c));fprintf('最大偏差： %f\n',d);d=sqrt(c'*c);fprintf('均方误差： %f',d);b=linspace(min(x),max(x),500);c=polyval(a,b);plot(x,y,'g*',b,c,'r');end结果：一次拟合多项式最大偏差： 2.700000 均方误差： 3.271085 拟合多项式曲线：y=6.5500x-12.5000二次拟合多项式最大偏差： 1.950000 均方误差： 2.906888 拟合多项式曲线：y=0.1875x2+4.6750x-8.7500四、解P89页第1题在Matlab中的截图如下：y=6.5500x-12.5000y=0.1875x2+4.6750x-8.7500其中a 为拟合多项式的系数阵列。

通过最小二乘法拟合多项式得到的多项式接近与实验值（代入后）。

作业3：用龙贝格算法计算,dx x x ⎰1)sin(使截断误差不超过0.5e-6。

目的和意义；在一元函数的积分学中,我们已经熟知,若函数f(x)在区间［a, b ］上连续且其原函数为F(x) ,则可用牛顿―莱布尼兹公⎰-=ba b F a F dx x f )()()(来求定积分。

牛顿―莱布尼兹公式虽然在理论上或在解决实际问题中都起了很大的作用,但它并不能完全解决定积分的计算问题。

当被积函数f(x)在区间［a,b ］上连续时,要使得复合梯形公式比较精确地代替定积分⎰-=bab F a F dx x f )()()( 可将分点(即基点)加密,也就是将区间［a,b ］细分,然后利用复合梯形公式求积将收敛缓慢的梯形值序列加工成收敛迅速的n n n R C S ,,，这种加速方法称为龙贝格算法。

它可以把复杂的积分比较精确的求出来。

一、龙贝格算法Matlab 程序如下：function z=romberg(ff,a,b,epsi,k) %用龙贝格算法计算积分 if nargin==4 k=16;end %如果只有4个传入参数，则将二分次数k 设置为默认值16 t=zeros(1,k+1); format long;function y=f(x) %嵌套函数，用来判断积分区域内的奇点y=ff(x); %若ff 参数为@(x)sin(x)/x,则有y=sin(x)/x if isnan(y)==true||isinf(y)==truedisp('请确保所给函数在积分区域内无奇点'); end endt(1)=(b-a)/2*(f(a)+f(b)); for l=1:k sum=0;for i=1:(2^(l-1))sum=sum+f(a+((2*i-1)*(b-a)/(2^l))); endt(l+1)=t(l)/2+(b-a)/(2^l)*sum; end %计算梯形序列 for m=1:k temp=t;for i=m+1:k+1 %将新序列放入t(m+1)--t(k+1)中 t(i)=((4^m)*temp(i)-temp(i-1))/(4^m-1); endif abs(t(m+1)-t(m))<epsi break; end endz=t(m+1);fprintf('截断误差为 %d\n',t(m+1)-t(m)); end二、用龙贝格算法计算P120页第5题：计算dx x x1)sin(,使截断误差不超过0.5e-6。

结果：截断误差为 6.632355e-008 =0.946083070387222三、程序流程图如下三、Matlab中的截图如下：二分次数缺省时：指定二分次数为10次时：积分区域内有奇点时：使用Matlab自带函数计算积分时：。

合集下载

计算方法课程上机实验报告2

计算方法课程上机实验报告2

计算方法课程上机实验二姓名 彭鑫 学号 2704302012 得分实验地点: 科研楼机房 学时: 指导教师: 李浩 时间:1. 实验目的掌握高斯消去法的基本思路和迭代步骤；培养编程与上机调试能力。

2. 实验内容利用列主元高斯消去法解方程组1231231230.101 2.304 3.555 1.1831.347 3.712 4.623 2.1372.835 1.072 5.643 3.035x x x x x x x x x ++=⎧⎪-++=⎨⎪-++=⎩要求给出相关理论基础和算法描述。

3. 算法介绍或方法基础高斯消去法对于一般方程组，通过同解变换将其化为上三角方程组这种特形式，然后进行求解.4. 程序function X=backsub(A,B)n=length(B);X=zeros(n,1);X(n)=B(n)/A(n,n);for k=n-1:-1:1X(k)=(B(k)-A(k,k+1:n)*X(k+1:n))/A(k,k);endfunction X=uptrbk(A,B)[N N]=size(A);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0breakendfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendAug;X=backsub(Aug(1:N,1:N),Aug(:,N+1));5. 实验结果X=uptrbk(A,B)X =-0.39820.01380.33516. 结果分析与解释从以上的结果可以知道1230.3982,0.0138,0.3351x x x =-==，高斯消去法计算量小，速度快，相比其他常规方法更有优势。

计算方法上机实验报告-MATLAB

计算方法上机实验报告-MATLAB

《计算方法》实验报告指导教师：学院：班级：团队成员：一、题目例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解k x ，并使其满足8110k k x x ---<原理：在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()()00x f x ，引切线()()()1000:'l y f x f x x x =+-其与x 轴相交于点：()()0100 'f x x x f x =-，进一步,过曲线()y f x =的点()()11x f x , 引切线()()()2111: 'l y f x f x x x =+-其与x 轴相交于点：()()1211 'f x x x f x =-如此循环往复，可得一列逼近方程()0f x =精确解*x 的点01k x x x ,,,,,其一般表达式为：()()111 'k k k k f x x x f x ---=-该公式所表述的求解方法称为Newton 迭代法或切线法。

程序：function y=f(x)%定义原函数y=x^3-x-1;endfunction y1=f1(x0)%求导函数在x0点的值syms x;t=diff(f(x),x);y1=subs(t,x,x0);endfunction newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数while abs(x1-x0)>=tolx0=x1;x1=x0-f(x0)/f1(x0);k=k+1;endfprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1);fprintf('迭代次数为k=%d\n',k);end结果：二、题目例3.7试利用Jacobi 迭代公式求解方程组1234451111101112115181111034x x x x ----⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪---⎝⎭⎝⎭⎝⎭ 要求数值解()k X 满足()4210k --≤X X ，其中T (1,2,3,4)=X 为方程组的精确解。

(完整word版)计算方法上机报告-备份

(完整word版)计算方法上机报告-备份

计算方法上机作业报告姓名：李小盼班级：计算方法B3班学号：6230489477专业：热能工程2016年《计算方法B》上机题目一.计算机语言要求使用语言不做限制，可以使用C、C++、FORTRAN、VC、VB、C#、Matlab、PHP、JavaScript等语言。

二.实习报告内容上机题目完成后，必须交一份上机报告。

上机报告中应对每道题目包括：（1）上机题目内容;（2）详细说明实现的思想、算法依据、算法实现的结构;（3）详细完整的源程序，并附相关的注释说明;（4）给出必要的计算结果（若数据量大，则只需列出主要的数据内容），并对结果进行分析;（5）对上机中出现的问题进行分析总结；三.实习报告要求1.提供一份完整的上机报告的电子文档；然后再提供一份与电子文档对应的纸质上机报告。

2.电子文档中提供上机过程中的所有源程序、输出数据，以及可以运行的文件。

如果程序的运行环境特殊，请附上运行程序所需要的软件环境。

3.上机报告严禁抄袭，如发现有抄袭现象，所有涉及抄袭的上机报告将被以作弊处理，并按零分处理,不再另行通知。

4.上机报告电子版一、二、三班发送到邮箱:xjtujsff@,上机作业纸面作业请送到:理科楼338。

上机实习题目1.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。

在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。

已探测到一组等分点位置的深度数据(单位:米)如下表所示：(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；解：1.1 实现思想选用曲线拟合数据点时，一方面要满足插值条件，即保证所得曲线经过以上所有数据点，另一方面也要保证曲线具有足够的“光滑性”，故这里采用三次样条插值法构造插值函数。

为构成封闭方程组，边界条件选取自然三次样条，以获得三次样条插值的关键参数M0和M n。

1.2 算法的依据以及结构三次样条插值具有较好的光滑性，可以用于对数据点的光滑拟合。

(完整word版)计算方法A上机实验报告

(完整word版)计算方法A上机实验报告

计算方法A上机实验报告姓名：苏福班级：硕4020 学号：3114161019一、上机练习目的1）复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

2）利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。

二、上机练习任务1）利用计算机语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。

2）掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。

3）写出上机练习报告。

三、上机题目1. 共轭梯度法求解线性方程组。

（第三章）2. 三次样条插值（第四章）3. 龙贝格积分（第六章）4. 四阶龙格－库塔法求解常微分方程的初值问题四、上机报告题目1：共轭梯度法求解线性方程组1．算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的问题。

从任意给定的初始点出发，沿一组关于矩阵A共轭的方向进行线性搜索，在无舍入误差的假定下，最多迭代n 次（其中n 为矩阵A 的阶数），就可求得二次函数的极小值，也就求得了线性方程组Ax b =的解。

定理：设A 是n 阶对称正定矩阵，则x *是方程组Ax b =的解得充分必要条件是x *是二次函数1()2TT f x x Ax b x =-的极小点，即 ()()min nx R Ax b f x f x **∈=⇔=共轭梯度法的计算公式：(0)(0)(0)()()()()(1)()()(1)(1)(1)()()()(1)(1)()k T k k k T k k k k k k k k T k k k T k k k k k d r b Ax r d d Ad xx d r b Ax r Ad d Ad d r d ααββ++++++⎧==-⎪⎪=⎪⎪=+⎪⎨=-⎪⎪⎪=-⎪⎪=+⎩2. 程序框图（1）编写共轭梯度法求解对称正定矩阵的线性方程组见附录（myge.m）：function x=myge(A,b)输入对称正定矩阵及对应的列向量，初始向量设为0，精度取为810 。