计算方法上机题问题详解

计算方法上机题答案

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. (25分)计算积分
dx x x I n n ⎰
+=105， n=0,1,2,…,20 若用下列两种算法 (A) n I I n n 151+
-=- (B) ⎪⎭
⎫ ⎝⎛-=-n n I n I 1511 试依据积分I n 的性质及数值结果说明何种算法更合理。

2. (25分)求解方程f(x)=0有如下牛顿迭代公式
)()(111---'-
=n n n n x f x f x x ， n ≥1，x 0给定 (1) 编制上述算法的通用程序，并以ε<--1n n x x （ε为预定精度）作为终止迭
代的准则。

(2)
利用上述算法求解方程 0cos :)(=-=x x x f
这里给定x 0=π/4，且预定精度ε=10-10。

3. (25分) 已知)3()(x x e x e x f -=，
(1)
利用插值节点x 0=1.00，x 1=1.02，x 2=1.04，x 3=1.06，构造三次Lagrange 插值公式，由此计算f(1.03)的近似值，并给出其实际误差； (2) 利用插值节点x 0=1，x 1=1.05构造三次Hermite 插值公式，由此计算f(1.03)的近
似值，并给出其实际误差。

4. (25分) 利用Romberg 算法计算积分
⎰
+4802cos 1dx x
精确到10
－4
总体要求：打印各题的程序代码及数值结果。

计算方法B上机报告

西安交通大学计算方法B上机报告班级：XXXXXXXX姓名：XXX学号：XXXXXXXX2015/12/13 Sunday目录1题目一 (1)1.数值计算 (1)2.实现思想 (1)3.源程序 (1)4.计算结果 (2)5.分析总结 (2)2题目二 (3)1.数据近似 (3)2.实现思想 (3)3.源程序 (5)4.计算结果 (6)5.分析总结 (7)3题目三 (8)1.数据拟合 (8)2.实现思想 (8)3.源程序 (9)4.计算结果 (12)5.分析总结 (14)4题目四 (15)1.非线性方程求解 (15)2.实现思想 (15)3.源程序 (16)4.计算结果 (17)5.分析总结 (17)5题目五 (18)1.线性方程组求解。

(18)2.文件格式说明 (18)3.实现思想 (19)4.源程序 (20)5.计算结果 (23)6.分析总结 (25)7.实际问题 (25)1 题目一1. 数值计算计算以下和式：0142118184858616nn S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求：(1)若保留11个有效数字，给出计算结果，并评价计算的算法；(2)若要保留30个有效数字，则又将如何进行计算。

2. 实现思想对于(1)可用迭代的方法进行处理。

通过do-while 循环使和式从0开始计算直到结果满足有效数字即可。

在循环中通过比较本次和上一次结果之差的绝对值与相应有效数字的大小的1/2，即可构成循环中止条件。

对于(2)可用同样的方法进行计算，然而由于所保留的有效数字较大，此时若不对上述算法进行改善的话，由于每一步所得计算结果都是其真实值，此时可能会存在有效数字丢失的情况，从而影响计算的准确性。

因此我们需要对程序中的计算精度进行控制。

在Maltlab 中可利用digits()函数对运算精度进行规定，为防止出现有效位的损失，在实际程序中需将精度提高几位。

而对迭代中每一步所得到的结果则利用vpa()函数进行限定，这样一来我们对每一个运算都控制了精度，而不只是单纯控制了结果的精度3. 源程序使用Matlab 所编写程序如下所示：(1)clear;clc;sd = 11; %所要保留有效数字 n = 0; %迭代次数S1 = 0; %当前n 时计算结果 S2 = 0; %n-1时计算结果 while 1eq = (4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6))/16^n; S1 = S1 + eq;if abs(S1-S2)<0.5*10^(-(sd-1)) %迭代是否终止判断条件 break endS2 = S1; n = n+1; endS=vpa(S1,sd)(2)clear;clc;sd = 30; %所要保留有效数字digits(sd+2);%控制精度n = 0; %迭代次数S1 = vpa(0); %当前n时计算结果S2 = vpa(0); %n-1时计算结果while 1eq = vpa((4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6))/16^n);S1 = vpa(S1 + eq);if abs(S1-S2)<vpa(0.5*10^(-(sd-1))) %迭代是否终止判断条件breakendS2 = S1;n = n+1;endS=vpa(S1,sd)4.计算结果(1)S =3.1415926536(2)S =3.141592653589793238435711679225.分析总结从计算结果可以看出与准确值3.1415926535897932384357116792180407…相比(1)和(2)所得结果分别是准确值保留11位和30位有效数字的结果，即利用本程序所得结果准确。

数值计算方法上机实习题答案

1．设⎰+=105dx x x I nn ，（1）由递推公式n I I n n 151+-=-，从0I 的几个近似值出发，计算20I ；解：易得：0I =ln6-ln5=0.1823,程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为：20I = -3.0666e+010（2）粗糙估计20I ，用n I I n n 515111+-=--，计算0I ；因为 0095.056 0079.010********≈<<≈⎰⎰dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为：I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endI0I = 0.0083（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。

首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。

并记nn n I I E '-=，则有01)5(5E E E n n n -==-=- 。

因为=20E 20020)5(I E >>-，所此递推式不可靠。

而在第二种递推式中n n E E E )51(5110-==-= ，误差在缩小，所以此递推式是可靠的。

出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。

2．求方程0210=-+x e x 的近似根，要求41105-+⨯<-k k x x ，并比较计算量。

（1）在[0，1]上用二分法；程序：a=0;b=1.0;while abs(b-a)>5*1e-4if exp(c)+10*c-2>0b=c;else a=c;endendc结果：c =0.0903（2）取初值00=x ，并用迭代1021x k e x -=+；程序：x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果：x =0.0905（3）加速迭代的结果；程序：x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果：x =0.0995（4）取初值00=x ，并用牛顿迭代法；程序：x=0;a=0;b=1;while abs(b-a)>5*1e-4x=x-(exp(x)+10*x-2)/(exp(x)+10);b=x;endx结果：x =0.0905（5）分析绝对误差。

浙大计算方法上机报告.

学号：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后运行下面的命令，直接输出结果。

计算方法上机题目

计算⽅法上机题⽬⽬录1.计算⽅法A 上机作业 (1)上机练习⽬的 (1)上机练习任务 (1)计算⽅法A 上机题⽬ (1)程序设计要求 (1)上机报告要求 (1)2.QR 分解法求解线性⽅程组 (2)计算原理 (2)程序框图 (7)计算实习 (8)Matlab代码 (8)3.共轭梯度法求解线性⽅程组 (10)计算原理 (10)程序框图 (11)计算实习 (12)Matlab代码 (12)4.三次样条插值 (14)计算原理 (14)程序框图 (16)计算实习 (17)Matlab代码 (17)5.四阶龙格－库塔法求解常微分⽅程的初值问题 (21)计算原理 (21)程序框图 (22)计算实习 (23)Matlab代码 (23)1.计算⽅法A 上机作业上机练习⽬的复习和巩固数值计算⽅法的基本数学模型，全⾯掌握运⽤计算机进⾏数值计算的具体过程及相关问题。

利⽤计算机语⾔独⽴编写、调试数值计算⽅法程序，培养学⽣利⽤计算机和所学理论知识分析解决实际问题的能⼒。

上机练习任务利⽤计算机语⾔编写并调试⼀系列数值⽅法计算通⽤程序，并能正确计算给定题⽬，掌握调试技能。

掌握⽂件使⽤编程技能，如⽂件的各类操作，数据格式设计、通⽤程序运⾏过程中⽂件输⼊输出运⾏⽅式设计等。

写出上机练习报告。

计算⽅法A 上机题⽬1. QR 分解⽅法求解线性⽅程组。

（第⼆章）2. 共轭梯度法求解线性⽅程组。

（第三章）3. 三次样条插值（第四章）4. 四阶龙格－库塔法求解常微分⽅程的初值问题程序设计要求1. 程序要求是通⽤的，在程序设计时要充分考虑哪些变量应该可变的。

2. 程序要求调试通过。

上机报告要求报告内容包括：●每种⽅法的算法原理及程序框图。

●程序使⽤说明。

●算例计算结果。

2. QR 分解法求解线性⽅程组计算原理当nx R∈是任意给定的⾮零向量，nv R∈是任意给定的单位向量，则存在初等反射阵2THI u u=-,使得H xvσ=，其中σ为常数，当取单位向量2x v ux vσσ-=-时，由u 确定的矩阵H 必定满⾜H xvσ=，所以在计算过程中取u 的值为上述值。

计算方法上机作业

计算方法上机报告姓名：学号：班级：上课班级：说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。

1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：①若只需保留11个有效数字，该如何进行计算； ②若要保留30个有效数字，则又将如何进行计算；（1）算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：142111416818485861681n n na n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭；2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1.;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n ；2.for 0,1,2,,n i =⋅⋅⋅if 10mt -≤end;3.for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。

计算方法上机题目

目录1.计算方法A 上机作业 (1)上机练习目的 (1)上机练习任务 (1)计算方法A 上机题目 (1)程序设计要求 (1)上机报告要求 (1)2.QR 分解法求解线性方程组 (2)计算原理 (2)程序框图 (8)计算实习 (9)Matlab代码 (11)3.共轭梯度法求解线性方程组 (13)计算原理 (13)程序框图 (14)计算实习 (15)Matlab代码 (19)4.三次样条插值 (21)计算原理 (21)程序框图 (23)计算实习 (24)Matlab代码 (26)5.四阶龙格－库塔法求解常微分方程的初值问题 (30)计算原理 (30)程序框图 (31)计算实习 (32)Matlab代码 (33)1.计算方法A 上机作业上机练习目的❑ 复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

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

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

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

•写出上机练习报告。

计算方法A 上机题目1. QR 分解方法求解线性方程组。

（第二章）2. 共轭梯度法求解线性方程组。

（第三章）3. 三次样条插值（第四章）4. 四阶龙格－库塔法求解常微分方程的初值问题程序设计要求1. 程序要通用的，在程序设计时要充分考虑哪些变量应该可变的。

2. 程序要求调试通过。

上机报告要求报告容包括：● 每种方法的算法原理及程序框图。

● 程序使用说明。

● 算例计算结果。

2. QR 分解法求解线性方程组计算原理当n x R ∈是任意给定的非零向量，n v R ∈是任意给定的单位向量，则存在初等反射阵2T H I uu =-,使得Hx v σ=，其中σ为常数，当取单位向量2x vu x vσσ-=-时，由u 确定的矩阵H 必定满足Hx v σ=，所以在计算过程中取u 的值为上述值。