当前位置：文档之家 › 西安交通大学计算方法上机作业

西安交通大学计算方法上机作业

  • 格式：doc
  • 大小：258.50 KB
  • 文档页数：19

计算方法上机作业1.对以下和式计算：0142118184858616n n S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若只需保留11个有效数字，该如何进行计算； (2)若要保留30个有效数字，则又将如何进行计算；（1）解题思想和算法实现：根据保留有效位数的要求，可以由公式得出计算精度要求。

只需要很少内存，时间复杂度和d 呈线性，不需要高浮点支持。

先根据while 语句求出符合精度要求的n 值的大小，然后利用for 语句对这n 项进行求和，输出计算结果及n 值大小即可。

（2）matlab 源程序：保留11位有效数字时； clear clcformat long n=0;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); while sum>=5*10^(-11); n=n+1;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); endfor i=0:n-1;sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); endvpa(sum,11) n保留30位有效数字时； clear clcformat long n=0;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); while sum>=5*10^(-30); n=n+1;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); endfor i=0:n-1;sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); endvpa(sum,30) n（3）实验结果分析图1.1 保留11位有效数字的n值及计算结果图图1.2 保留30位有效数字的n值及计算结果图由计算结果可知，通过合理的误差控制，分别通过7次和22次循环，可以实现题目所要求的精确度。

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

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

已探测分点0 1 2 3 4 5 6深度9.01 8.96 7.96 7.97 8.02 9.05 10.13分点7 8 9 10 11 12 13深度11.18 12.26 13.28 13.32 12.61 11.29 10.22分点14 0深度9.15 7.90 7.95 8.86 9.81 10.80 10.93(1)请用合适的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；（1）解题思想和算法原理给定区间[a, b]一个分划⊿：a=x0<x1<…<x N=b若函数S(x)满足下列条件：1)S(x)在每个区间[x i, x j]上是不高于3次的多项式。

2)S(x)及其2阶导数在[a, b]上连续。

则称S(x)使关于分划⊿的三次样条函数。

（2）matlab源程序：clc,clearx=0:1:20;y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.9 7.95 8.86 9.81 10.80 10.93];n=length(x);l(1)=0;m(n)=0;h=diff(x);df=diff(y)./diff(x);d(1)=0;d(n)=0;for j=2:n-1l(j)=h(j)/(h(j-1)+h(j));m(j)=h(j-1)/(h(j-1)+h(j));d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j));endm=m(2:end);u=diag(m,-1);r=diag(l,1);a=diag(2*ones(1,n));A=u+r+a;M=inv(A)*d';syms gfor j=1:n-1s(j)=M(j)*(x(j+1)-g)^3/(6*h(j))+M(j+1)*((g-x(j))^3/(6*h(j)))+(y(j)-M(j)*h(j)^2/6)*(x(j+1) -g)/h(j)+(y(j+1)-M(j+1)*h(j)^2/6)*(g-x(j))/h(j);endsr=0;for j=1:n-1df=diff(s(j),g);warning off all;q=int(sqrt(1+df.^2),g,j-1,j);r=r+q;endL=vpa(r,8);disp('the length of the label is L=');disp(L);for j=1:n-1S(j,:)=sym2poly(s(j));endfor j=1:n-1x1=x(j):0.1:x(j+1);y1=polyval(S(j,:),x1);if j==1y2=y1;elsefor i=1:11k=(j-1)*10+i;y2(k)=y1(i);endendendx2=x(1):0.1:x(n);plot(x,y,'o')gridhold onplot(x2,y2,'r')（3）实验结果分析图2.1 铺设河底电缆长度图2.2 铺设河底光缆的曲线图由三次样条插值得出的函数曲线的长度和即铺设河底电缆的长度为26.498514。

为了提高插值精度，用三次样条插值可以增加插值节点的办法来满足要求，而且在给定节点数的条件下，三次样条插值的精度要优于多项式插值以及线性分段插值，虽然舍弃了降低误差这个优点，但是其曲线的光滑性要好一些。

3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这（1）解题思想和数学原理：对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。

又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。

最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。

对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为∑==mj j j x a x y 0)()(ϕ特别的，取)(x j ϕ为多项式jj x x =)(ϕ (j=0, 1,…，m )则根据最小二乘法原理，可以构造泛函∑∑==-=ni mj i j j i m x a f a a a H 110))((),,,(ϕ令0=∂∂ka H(k=0, 1,…，m )则可以得到法方程⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ求该解方程组，则可以得到解ma a a ,,,10 ，因此可得到数据的最小二乘解∑=≈mj j j x a x f 0)()(ϕ（2）matlab 源程序：x=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]; %给出题目数据（时间）y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16]; %给出题目数据（温度）plot(x,y, 'm*') %画出各个离散数据点 hold onfor n=2:4; %2、3、4代表拟合函数最高项次数 alltemp=25; % alltemp 代表数据点总共有25个A=zeros(n+1,n+1); %定义初始正规方程组的系数矩阵A C=ones(n+1,1); %定义初始正规方程组的系数向量C D=zeros(n+1,1); %定义初始正规方程组的向量D for i=1:n+1for j=1:n+1for k1=1: alltempA(i,j)=A(i,j)+(x(k1).^(i-1+j-1)); end endfor k2=1: alltempD(i,1)=D(i,1)+(x(k2).^(i-1)).*(y(k2)); endend %以上为计算出正规方程组矩阵A 、D 的所有元素的程序 tol=1.0e-12; maxit=1000;C=bicg(A,D,tol,maxit); %使用bicg 迭代算出正规方程组的系数向量C p=0; %误差分量 E=0; %误差总量 if n==2b=0:24;f=C(1)+C(2).*b+C(3).*(b.^2); p=y(b+1)-f; for v=1:25E=E+(p(v)).^2; endplot(b,f, 'r-')end %以上是对2阶拟合函数的图形处理与误差计算 if n==3b=0:24;f=C(1)+C(2).*b+C(3).*(b.^2)+C(4).*(b.^3);p=y(b+1)-f;for v=1:25E=E+(p(v)).^2;endplot(b,f, 'g-')end %以上是对3阶拟合函数的图形处理与误差计算if n==4b=0:24;f=C(1)+C(2).*b+C(3).*(b.^2)+C(4).*(b.^3)+C(5).*(b.^4);p=y(b+1)-f;for v=1:25E=E+(p(v)).^2;endplot(b,f, 'b-')end %以上是对4阶拟合函数的图形处理与误差计算C,Eendn=2; %重新对n赋值，进行指数函数拟合A=zeros(n+1,n+1); %重新对A矩阵赋初值C=zeros(n+1,1); %重新对C向量赋初值D=zeros(n+1,1); %重新对D向量赋初值for i=1:n+1for j=1:n+1for k=1: alltempA(i,j)=A(i,j)+(x(k).^(i-1+j-1));endendfor l=1: alltempD(i,1)=D(i,1)+((x(l).^(i-1)).*(log(y(l))));endend %计算出A矩阵、D向量各元素数值C=bicg(A,D,tol,maxit); %利用bicg迭代求解系数b=0:24;p=0;E=0;f=exp(C(1)+C(2).*b+C(3).*(b.^2));p=y(b+1)-f;for v=1:25E=E+(p(v)).^2;endplot(b,f, 'c-') %对指数拟合函数进行图形处理和误差计算b=-C(3);c=C(2)/(2*b);a=exp(b*(c^2)+C(1)); %算出题设要求的指数拟合函数的各个系数a,b,c,Egrid onlegend('测量数据','二次函数','三次函数','四次函数','指数拟合','Location','NorthWest') hold off %hold on与hold off联合使用，表示将各个曲线画在同一个图中图3.1 二次曲线拟合系数与2范数误差图3.2 三次曲线拟合系数与2范数误差图3.3 四次曲线拟合系数与2范数误差图3.4 指数曲线拟合系数与2范数误差图3.5 数据原始点与拟合曲线对比图（3）实验结果分析：根据国家有关规定，用的是02、08、14、20时四个观测时次的数据做平均，最有代表性。

合集下载

计算方法上机报告_董晓壮

计算方法上机报告_董晓壮

计算方法A上机报告学院（系）：电气工程学院学生姓名：陶然学号：授课老师：完成日期：2019年12月03日西安交通大学Xi'an Jiaotong University目录1 QR分解法求解线性方程组 (2)1.1 算法原理 (2)1.1.1 基于吉文斯变换的QR分解 (2)1.1.2 基于豪斯霍尔德变换的QR分解 (3)1.2 程序流程图 (4)1.2.1 基于吉文斯变换的QR分解流程图 (4)1.2.2 基于豪斯霍尔德变换的QR分解流程图 (5)1.3 程序使用说明 (5)1.3.1 基于吉文斯变换的QR分解程序说明 (5)1.3.2 基于豪斯霍尔德变换的QR分解程序说明 (7)1.4 算例计算结果 (8)2 共轭梯度法求解线性方程组 (10)2.1 算法原理 (10)2.2 程序流程图 (10)2.3 程序使用说明 (11)2.4 算例计算结果 (12)3 三次样条插值 (14)3.1 算法原理 (14)3.2 程序流程图 (16)3.3 程序使用说明 (17)3.4 算例计算结果 (19)4 龙贝格积分 (21)4.1 算法原理 (21)4.2 程序流程图 (22)4.3 程序使用说明 (23)4.4 算例计算结果 (24)结论 (26)1 QR 分解法求解线性方程组1.1 算法原理矩阵的QR 分解是指，可以将矩阵A 分解为一个正交矩阵Q 和一个上三角矩阵R 的乘积，实际中，QR 分解经常被用来解决线性最小二乘问题，分解情况如图1.1所示。

=⨯图1.1 QR 分解示意图本次上机学习主要进行了两个最基本的正交变换—吉文斯（Givens ）变换和豪斯霍尔德（Householder ）变换，并由此导出矩阵的QR 分解以及求解线性方程组的的方法。

1.1.1 基于吉文斯变换的QR 分解吉文斯矩阵也称初等旋转阵，如式（1.1）所示，它把n 阶单位矩阵I 的第,i j 行的对角元改为c ，将第i 行第j 列的元素改为s ，第j 行第i 列的元素改为s −后形成的矩阵。

(完整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 。

计算方法上机题答案

计算方法上机题答案

2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量（1）二分法（局部，大图不太看得清，故后面两小题都用局部截图）（2）迭代法（3）牛顿法顺序消元法#include<stdio.h>#include<stdlib.h>#include<math.h>int main(){ int N=4,i,j,p,q,k; double m;double a[4][5]; double x1,x2,x3,x4; for (i=0;i<N ;i++ )for (j=0;j<N+1; j++ )scanf("%lf",&a[i][j]);for(p=0;p<N-1;p++) {for(k=p+1;k<N;k++){m=a[k][p]/a[p][p];for(q=p;q<N+1;q++)a[k][q]=a[k][q]-m*a[p][q];}}x4=a[3][4]/a[3][3];x3=(a[2][4]-x4*a[2][3])/a[2][2];x2=(a[1][4]-x4*a[1][3]-x3*a[1][2])/a[1][1];x1=(a[0][4]-x4*a[0][3]-x3*a[0][2]-x2*a[0][1])/a[0][0];printf("%f,%f,%f,%f",x1,x2,x3,x4);scanf("%lf",&a[i][j]); （这一步只是为了看到运行的结果）}运行结果列主元消元法function[x,det,flag]=Gauss(A,b)[n,m]=size(A);nb=length(b);flag='OK';det=1;x=zeros(n,1);for k=1:n-1 max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endendif max1<1e-10flag='failure';return;endif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det; endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n)if abs(A(n,n))<1e-10flag='failure';return;endx(n)=b(n)/A(n,n);for k=n-1:-1：1for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end运行结果：雅可比迭代法function y=jacobi(a,b,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>1e-4 x0=y;y=B*x0+f;n=n+1;endyn高斯赛德尔迭代法function y=seidel(a,b,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>10^(-4) x0=y;y=G*x0+f;n=n+1;endynSOR迭代法function y=sor(a,b,w,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);lw=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w;y=lw*x0+f;n=1;while norm(y-x0)>10^(-4) x0=y;y=lw*x0+f;n=n+1;endyn1.分段线性插值：function y=fdxx(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for i=1:mz=x(i);for j=1:n-1if z<x0(j+1)break;endendy(i)= y0(j)*(z-x0(j+1))/(x0(j)-x0(j+1))+y0(j+1)*(z-x0(j))/(x0(j+1)-x0(j));fprintf('y(%d)=%f\nx1=%.3fy1=%.3f\nx2=%.3fy2=%.3f\n\n',i,y(i),x0(j),y0(j),x0(j+1),y0(j+1));endend结果0.39404 0.38007 0.356932.分段二次插值：function y=fdec(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for i=1:mz=x(i);for j=1:n-1if z<x0(j+1)break;endendif j<2j=j+1;elseif (j<n-1)if (abs(x0(j)-z)>abs(x0(j+1)-z))j=j+1;elseif ((abs(x0(j)-z)==abs(x0(j+1)-z))&&(abs(x0(j-1)-z)>abs(x0(j+2)-z)))j=j+1;endendans=0.0;for t=j-1:j+1a=1.0;for k=j-1:j+1if t~=ka=a*(z-x0(k))/(x0(t)-x0(k));endendans=ans+a*y0(t);endy(i)=ans;fprintf('y(%d)=%f\n x1=%.3f y1=%.3f\n x2=%.3f y2=%.3f\n x3=%.3f y3=%.3f\n\n',i,y(i),x0(j-1),y0(j-1),x0(j),y0(j),x0(j+1),y0(j+1));endend结果为0.39447 0.38022 0.357253.拉格朗日全区间插值function y=lglr(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for t=1:mans=0.0;z=x(t);for k=1:np=1.0;for q=1:nif q~=kp=p*(z-x0(q))/(x0(k)-x0(q));endendans=ans+p*y0(k);endy(t)=ans;fprintf('y(%d)=%f\n',t,y(t));endend结果为0.39447 0.38022 0.35722function [p,S,mu] = polyfit(x,y,n)if ~isequal(size(x),size(y))errorendx = x(:);y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);endV(:,n+1) = ones(length(x),1,class(x));for j = n:-1:1V(:,j) = x.*V(:,j+1);end[Q,R] = qr(V,0);ws = warning;p = R\(Q'*y); warning(ws);if size(R,2) > size(R,1)warning;elseif warnIfLargeConditionNumber(R)if nargout > 2warning;elsewarning;endendif nargout > 1r = y - V*p;S.R = R;S.df = max(0,length(y) - (n+1));S.normr = norm(r);endp = p.';function flag = warnIfLargeConditionNumber(R) if isa(R, 'double')flag = (condest(R) > 1e+10);elseflag = (condest(R) > 1e+05);endx=[1,3,4,5,6,7,8,9,10];y=[10,5,4,2,1,1,2,3,4];p=polyfit(x,y,2);y=poly2sym(p);a=-p(1)/p(2)*0.5;xa=-p(2)/(2*p(1));min=(4*p(1)*p(3)-p(2)^2)/(4*p(1)); yxamin运行截图运行结果function [I] = CombineTraprl(f,a,b,eps)if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym( f),findsym(sym(f)),x1));endendI=I2;程序：>> [q]=CombineTraprl('(1-exp(-x))^0.5/x'10^-12,1)结果：q =1.8521欧拉方法function [x,y]=euler(fun,x0,xfinal,y0,n);if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i));end程序：[x,y]=euler('doty',0,1,1,10)结果：改进欧拉方法function [x,y]=eulerpro(fun,x0,xfinal,y0,n); if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;x(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end程序：>> [x,y]=eulerpro('doty',0,1,1,10)结果：经典RK法function [x,y]=RungKutta4(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; end程序：dyfun=inline('(2*x*y^-2)/3');[x,y]=RungKutta4(dyfun,[0,1],1,0.1)标准函数x(1)=0;h=0.1;for i=1:10y(i)=(1+x(i)^2)^(1/3); endxY结果：。

计算方法大作业1 克服Runge现象

计算方法大作业1 克服Runge现象


x3
x2
x
1
S1 ( x)
-0.34685
0.2086
0.073964
0.038462
S2 (x)
S (xi 0 ) S x(i 0 )

S
'
(xi

0) S
xi' (
0 )i

S
'
'
x(i

0)S
xi' ' (
0)
1 ,n2, . . . , 1
（1）
这里共有了 3n-3 个条件，再加上条件（2）中的 n+1 个插值条件，共有 4n-2 个条件，
因此还需要 2 个方程才能确定 S (x) .通常可在区间[a, b]的端点 a x0,b xn 上各加一个边

dn1

1
2


Mn


dn

（6）
2 1


2
2
2
1 M1 d1

M2


d2




n 1
2
n
1


M
n
1

dn1
n
n 2 M n dn
由式（1）内点拼接条件，可得
i M i1 2M i i M i1 d j i 1, 2,..., n 1
（3） （4）
其中
i

hi 1 hi1
, hi

i

hi hi 1