北京科技大学计算方法大作业
计算方法大作业
机械电子工程系 老师：廖福成
注：本文本只有程序题，证明题全部在手写已交到理化楼204了。
2. 证明方程 3
10x x --=
在[1,2]上有一实根*x ，并用二分法求这个根。要求31||10k k x x -+-<。请给出程
序和运行结果。
证明：
设f(x)=x3-x-1
则f(1)= -1，f(2)= 5，f(1)*f(2)= -5<0 因此，方程在[1,2]上必有一实根。
二分法求解程序：
%预先定义homework2.m文件如下：function lc=homework2(x)
lc=x^3-x-1;
在MALAB窗口运行：
clear
a=1;b=2;tol=10^(-3);N=10000;
k=0; fa=homework2(a); % f 需事先定义for k=1:N
p=(a+b)/2;
fp=homework2(p);
if( fp==0 || (b-a)/2 break end if fa*fp<0 b=p; else a=p; end end k,p 程序运行结果： k = 10 p = 1.325195312500000 3. 用Newton 迭代法求方程 32 210200x x x ++-= 的一个正根，计算结果精确到7位有效数字. 要求给出程序和运行结果. 解： 取迭代初值01x = ,并设32()21020f x x x x =++-,则 '2 ()3410f x x x =++． 牛顿迭代函数为 32'2 ()21020 ()()3410f x x x x x x x f x x x ?++-=-=-++ 牛顿迭格式为：3212 21020 3410k k k k k k k x x x x x x x +++-=-++ Matlab 程序如下： %定义zuoye3.m 文件 function x=zuoye3(fname,dfname,x0,e,N) if nargin<5,N=500;end if nargin<4,e=1e-7;end x=x0;x0=x+2*e;k=0; while abs(x0-x)>e&k x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x) end if k==N,warning('已达上限次数');end 在Matlab 窗口中执行： zuoye3(inline('x^3+2*x^2+10*x-20'),inline('3*x^2+4*x+10'),1,1e-7) 结果如下： 1.41176470588235 1.36933647058824 1.36880818861753 1.36880810782137 ans = 1.36880810782137 4. 用牛顿迭代法求方程310x x --=在01 x =附近的根. 要求给出程序和运行结 果. 解：令：3()1f x x x =--,则 '2 ()31f x x =-． 牛顿迭代函数为 3'2()1 ()()31f x x x x x x f x x ?--=-=- - 牛顿迭格式为： 312 1 31k k k k k x x x x x +--=-- Matalb 程序如下： %定义zuoye4.m 文件 function x=zuoye4(fname,dfname,x0,e,N) if nargin<5,N=500;end if nargin<4,e=1e-7;end x=x0;x0=x+2*e;k=0; while abs(x0-x)>e&k k=k+1; x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x) end if k==N,warning('已达上限次数');end 在Matlab窗口执行： zuoye4(inline('x^3-x-1'),inline('3*x^2-1'),1,1e-7) 结果如下： 1.50000000000000 1.34782608695652 1.32520039895091 1.32471817399905 1.32471795724479 1.32471795724475 ans = 1.32471795724475 6. 编写用全主元Gauss消去法解线性方程组的程序，并求解 12345123451234512345 123450.024******* 42334332416 34418x x x x x x x x x x x x x x x x x x x x x x x x x -+-+=??-++++=?? +++-=??-++++=??+-++=? 解： Matlab 程序如下： A=[2 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4] b=[11 14 4 16 18] function x=zuoye6(A,b) [n,v]=size(b); D=[A,b;eye(n),zeros(n,v)],[s1,m]=size(D); for k=1:(n-1) s=abs(A(k,k));p=k;q=k; for i=k:n for j=k:n if abs(A(i,j))>s s=abs(A(i,j));p=i;q=j; end end end if p>k t=D(k,:); D(k,:)=D(p,:); D(p,:)=t; end if q>k t1=D(:,k); D(:,k)=D(:,q); D(:,q)=t1; end h=D(k+1:n,k)/D(k,k); D(k+1:n,k+1:m)=D(k+1:n,k+1:m)-h*D(k,k+1:m); D(k+1:n,k)=zeros(n-k,1); end for k=n:-1:1 D(k,k:m)=D(k,k:m)/D(k,k); for r=1:k-1 D(r,:)=D(r,:)-D(r,k)*D(k,:); end end P=D(n+1:2*n,1:n);Q=D(1:n,n+1:m);x=P*Q 在Matlab窗口中执行： A=[0.02 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4]; b=[11 14 4 16 18]'; zuoye6(A,b) 运行结果如下： x = 2.94117647058824 -3.82352941176472 1.00000000000000 0.94117647058824 5.94117647058824 7. 用追赶法解线性方程组 1 2 3 4 5 210001 121000 012100 001210 000120 x x x x x -?? ???? ?? ????--?? ???? ?? ???? = -- ?? ???? --?? ???? ?? ???? - ???? ?? 要求给出程序和运行结果. 解： 10000210 01310000100210 0022121002 4 1000010012103 30 012135 001000014 4000124600 100 00 5 5A LU -???????? ???? ---????????--??????--??????===--???? ?? ????--??--???? ??-?????? ???? -????? ?? ? 于是有求解Ax=b 即为求解{ Ly b Ux y ==，式中b=(1 0 0 0 0)T 据 123451000 0110001202 010003 3000104 0400 15y y y y y ??????-???? ????????????-?????? =????????????-??????????????-???? y=112131415???? ?????????????????????? 据210003010024 00103500014600005-?? ???? -????-???? ??-?? ?? ? ???? ?12345x x x x x ????????????????=112131415?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ? ??? x=5623121316?? ?????????????????????????? Matlab 程序如下： %定义zuoye7.m文件 function x=zuoye7(a,b,c,d) a1=[0;a]; n=length(b);q=zeros(n,1);p=zeros(n,1); %LU分解 q(1)=b(1);for k=2:n,p(k)=a1(k)/q(k-1); q(k)=b(k)-p(k)*c(k-1); end %解Ly=d y=zeros(n,1);y(1)=d(1);for k=2:n, y(k)=d(k)-p(k)*y(k-1);end %解Ux=y x=zeros(n,1); x(n)=y(n)/q(n);for k=n-1:-1:1,x(k)=(y(k)-c(k)*x(k+1))/q(k);end x 在Matlab窗口中执行： a=[-1 -1 -1 -1]'; b=[2 2 2 2 2]'; c=[-1 -1 -1 -1]'; d=[1 0 0 0 0]'; x=zuoye7(a,b,c,d) 运行结果如下： x = 0.83333333333333 0.66666666666667 0.50000000000000 0.33333333333333 0.16666666666667 9. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解程组（编写程序） 12312312310+3142-103-531014x x x x x x x x x +=????+=???? ++=?? 取初值=T X （0）（0,0,0），精确到小数后面四位。 解： （1） Jacobi 迭代法的Matlab 程序： % x0为初始向量，ep 为精度，N 为最大次数，x 是近似解向量 A=[10 3 1;2 -10 3;1 3 10]; b=[14 -5 14];n=length(b);N=500;ep=1e-6;x0=zero(n,1); n=length(b);x0=zeros(n,1);x=zeros(n,1);k=0 while K x(i)=(b(i)-A(I,[1:i-1,i+1:n])*x0(1:i-1,i+1:n))/A(i,i); end if norm(x-x0,inf) if k==N,Warning(‘已达到迭代次数上限’);end disp([‘k=’,num2str(k)]),x 运行结果如下： k=15 x= 1.000000327906423 0.999999801006905 1.000000327906432 （2）Gauss-Seidel迭代法的Matlab程序： % x0为初始向量，ep为精度，N为最大次数，x是近似解向量 Format long;clear; A=[10 3 1;2 -10 3;1 3 10]; b=[14 -5 14];n=length(b);N=500;ep=1e-6;x0=zero(n,1);P=inf; %以下是Guass-Seidal迭代法程序 D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);dD=det(D); if dD==0 disp(‘请注意：因为对角矩阵D奇异，所以此方程组无解。’) else disp(‘请注意：因为对角矩阵D非奇异，所以此方程组无解。’) iD=inv(D-L);B_GS=ID*U;d_GS=iD*b; end k=0; while k x=B_GS*x0+d_GS; if norm(x-x0,P) x0=x;k=k+1; end if k==N,warning(‘已达到迭代次数上限’)；end disp([‘k=’,num2str(k)]),x,%D,U,L,B_GS,d_GS,%jx=A\b, 运行结果如下： 请注意：因为对角矩阵D 非奇异，所以此方程组有解。 k=9 x=1.000000043651052 1.000000031224555 0.999999986267528 10．编写幂法程序求矩阵 1 10.5110.250.50.252A ????=?? ???? 按模最大的特征值1λ和对应的特征向量1x 。6(10)ε-= 解：幂法程序如下： % A 为n 阶方阵，u0 为初始向量， %epsilon 为上限，max1 为循环次数。lambda 返回按模最大的特征值，u 返回 对应的特征向量。 clear;format long;clc;max1=1000; epsilon=1e-6; A=[1 1 0.5;1 1 0.25;0.5 0.25 2]; u0=[1;1;1]; u=u0;v0=u0; lambda=0;k=0;err=1;R=[]; while ((k v=A*u; [m,j]=max(abs(v)); m=v(j); u=v/m; err=abs(lambda-m); lambda=m; k=k+1; Temp=[k;m;v;u];R=[R Temp]; end k, lambda, u, disp(R), xlswrite('d:\\1.xls', R,'sheet1', 'b1') ; %[D,Lmd]=eig(A) %D 之列为A 的特征向量，Lmd 之对角线上为A 的特征值 运行结果如下： k = 23 lambda = 2.536527202038736 u = 0.748222175139906 0.649662222855908 1.000000000000000 11. 编写用原点位移加速反幂法程序求矩阵 ?? ??? ?????= 3 1- 2-1- 4 3 2- 3 7 A 最接近于 1.9q =的特征值和相应的特征向量。 10 (10)ε-= 解： 反幂法程序如下： % A 为n 阶方阵，v0,u0 为初始向量， %epsilon 为上限，max1 为循环次数，q 为近似特征值。lambda 为提高精度的 特征值，v 给出对应的特征向量。 clear;format long;max1=100; epsilon=1e-10; A=[7 3 -2;3 4 -1;-2 -1 3]; u0=[1;1;1]; u=u0;v=u0;q=1.9; lambda=0;k=0;err=1;mu=0.5; n=length(u0); z=zeros(n,1);A=A-q*eye(n); %LU 分解 n=length(u0);UU=zeros(n,n);L=eye(n,n); UU(1,:)=A(1,:); L(2:n,1)=A(2:n,1)/UU(1,1); for s=2:n UU(s,s:n)=A(s,s:n)-L(s,1:s-1)*UU(1:s-1,s:n); L(s+1:n,s)=(A(s+1:n,s)-L(s+1:n,1:s-1)*UU(1:s-1,s))/UU(s,s); end while ((k [m,j]=max(abs(v)); m=v(j); u=v/m; %解下三角方程组Lz=u z(1)=u(1);for j=2:n, z(j)=u(j)-L(j,1:j-1)*z(1:j-1); end %解上三角方程组Uv=z v(n)=z(n)/UU(n,n); for j=n-1:-1:1, v(j)=(z(j)-UU(j,j+1:n)*v(j+1:n))/UU(j,j); end err=abs(1/m-1/mu); k=k+1; mu=m; end disp(['迭代次数k= ',num2str(k)]) disp('要求的矩阵A 的特征值lambda='), lambda= q+1/m disp('要求的矩阵A 的特征向量u'), u 运行结果如下： 迭代次数k= 17 要求的矩阵A 的特征值lambda= lambda = 2.000000000009091 要求的矩阵A 的特征向量u u = 0.999999999973922 -0.999999999956935 1.000000000000000 12. 已知插值点 (-2.00,17.00), (0.00,1.00), (1.00,2.00), (2.00,17.00), f. 求三次插值多项式，并计算(0.6) 解： 其Matlab程序如下： function yy=zuoye10(x,y,xx) m=length(x);n=length(y); if m~=n, error('向量x与y的长度必须一致');end for i=1:n t=ones(1,length(xx)); for j=1:n if j~=i t=t.*(xx-x(j))/(x(i)-x(j)); end end s=s+t*y(i); end yy=s; 在Matlab 窗口中执行： x=[-2.00 0.00 1.00 2.00]; y=[17.00 1.00 2.00 17.00]; xx=0.6; zuoye10(x,y,xx) 结果如下： ans = 0.25600000000000 13．编制Newton 插值法的通用Matlab 程序，并求(0.5)f 的近似值. 已知的数值如下表所示 00.20.40.60.8 ()0.19950.39650.58810.7721.09461 i i x f x Newton 插值法的通用Matlab 程序： x=[0 0.2 0.4 0.6 0.8]';y=[0.1995 0.3965 0.5881 0.7721 0.9461]';xx=0.5; m=length(x);n=length(y); if m~=n, error('向量x 与y 的长度必须一致');end for j=2:n ，for i=j:n newton_chazhi(i,j)=(newton_chazhi(i-1,j-1)-newton_chazhi(i,j-1))/(x(1+i-j)-x(i)); end end yy= newton_chazhi(1,1); t=1;for i=2:n ，t=t*(xx-x(i-1));yy=yy+ newton_chazhi(i,i)*t; end %计算误差近似值 err=newton_chazhi(n,n);for i=1:n,err=err*(xx-x(i));end,err yy 运行结果为： err =-2.343750000006213e-06 yy=0.681187500000000 20．编写解常微分方程00 (,)|x x y f x y y y ='=??? =??的四阶龙格库塔法通用程序. 解： 四阶龙格库塔法通用程序： function [x,y]=zuoye20(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=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 x=x'; y=y'; 21. 利用上题的程序求解初值问题: 22,[0,1.2] 3(0)1 dy xy x dx y -?=∈???=? 的数值解, 并将结果与准确解y =. 取0.4h =. 解： Matlab 程序： function [x,y]=zuoye21(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=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 x=x' y=y' 在Matlab窗口中执行： dyfun=inline('2*x*y^(-2)/3'); xspan=[0,1.2]; y0=1; h=0.4; nark18(dyfun,xspan,y0,h) 结果如下： x = 0.40000000000000 0.80000000000000 1.20000000000000 y =1.00000000000000 1.05075062243354 1.17933175912031 1.34631542500313 计算值与准确值比较如下： 22. 编写四阶龙格库塔法程序，并求解洛伦兹型系统（不许调用现成的龙格库塔法程序软件包）： ??? ??+-=+-=++-=xy bz z sxz qy rx y yz y x x μετσ 解：为了便于求解，这里首先为各个参数赋予初值，实际中根据需要修改程序中的参数即可。 取25.0=σ，06.0=τ，4.0=b ，5.0=ε，120=r ，3.1=q ，5.1=s ，20-=μ.取初始点为(0.005,0.4596,0.1146)-. 程序如下： %预定义homework22.m 文件 function f=homework22(t,y) 计算方法 姓名： 学号： 班级： 指导教师： 目录 作业1 (1) 作业2 (5) 作业3 (8) 作业4 (10) 作业5 (14) 作业6 (16) 作业7 (17) 作业1 1、分别用不动点迭代与Newton 法求解方程 -+=x 2x e 30的正根与负根。 解： (1)不动点迭代 a.原理： 将 230x x e -+=变型为1()k k x g x +=进行迭代，直到为止 变型后为有两种形式：和 b.程序：初值为1 形式: x=zeros(100,1); tol=1; i=1; x(1)=1; whiletol>=10e-6; disp(x(i)) x(i+1)=log(2*x(i)+3); tol=abs(x(i+1)-x(i)); i=i+1; end disp(i-1); 形式: x=zeros(100,1); tol=1; i=1; x(1)=1; whiletol>=10e-6; disp(x(i)) x(i+1)=(exp(x(i))-3)/2; tol=abs(x(i+1)-x(i)); i=i+1; end disp(i-1); c.运行结果： 初值为1 (23) 1ln k x k x ++=6 110k k x x -+-<13 2 k x k e x +-= (23)1ln k x k x ++=132 k x k e x +-= 迭代次数：11 迭代次数：9 (2)Nexton法 a.原理： 令 () () 1' k k k k f x x x f x + =-得到迭代公式为： () 1 23 2 k k x k k k x x e x x e + -+ =- - b.程序：初值为0 x=zeros(100,1); tol=1; i=1; x(1)=0; whiletol>=10e-6; disp(x(i)) x(i+1)=x(i)-((2*x(i)-exp(x(i))+3)/(2-exp(x(i)))); tol=abs(x(i+1)-x(i)); i=i+1; end disp(i-1); 初值为1 x=zeros(100,1); tol=1; i=1; x(1)=1; whiletol>=10e-6; disp(x(i)) x(i+1)=x(i)-((2*x(i)-exp(x(i))+3)/(2-exp(x(i)))); tol=abs(x(i+1)-x(i)); i=i+1; end disp(i-1) a=x(i-1); b=2*a-exp(a)+3; disp(b); c.运行结果： 初值为0 目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20) 4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28) 数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10） 1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10） 2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量 （1）二分法 （局部，大图不太看得清，故后面两小题都用局部截图） （2）迭代法 （3）牛顿法 顺序消元法 #include for(k=p+1;k 一 判断题（20分） （×）1. 节点的位置依赖于形态，而并不依赖于载荷的位置 （√）2. 对于高压电线的铁塔那样的框架结构的模型化处理使用梁单元 （×）3. 不能把梁单元、壳单元和实体单元混合在一起作成模型 （√）4. 四边形的平面单元尽可能作成接近正方形形状的单元 （×）5. 平面应变单元也好，平面应力单元也好，如果以单位厚来作模型化 处理的话会得到一样的答案 （×）6. 用有限元法不可以对运动的物体的结构进行静力分析 （√）7. 一般应力变化大的地方单元尺寸要划的小才好 （×）8. 所谓全约束只要将位移自由度约束住，而不必约束转动自由度 （√）9. 同一载荷作用下的结构，所给材料的弹性模量越大则变形值越小 （√）10一维变带宽存储通常比二维等带宽存储更节省存储量。 二、填空（20分） 1．平面应力问题与薄板弯曲问题的弹性体几何形状都是 薄板 ，但前者受力特点是： 平行于板面且沿厚度均布载荷作用 ，变形发生在板面内； 后者受力特点是： 垂直于板面 的力的作用，板将变成有弯有扭的曲面。 2．平面应力问题与平面应变问题都具有三个独立的应力分量： σx ，σy ，τxy ，三个独立的应变分量：εx ，εy ，γxy ，但对应的弹性体几何形状前者为 薄板 ，后者为 长柱体 。3．位移模式需反映 刚体位移 ，反映 常变形 ，满足 单元边界上位移连续 。 4．单元刚度矩阵的特点有：对称性 ， 奇异性 ，还可按节点分块。 5．轴对称问题单元形状为：三角形或四边形截面的空间环形单元 ，由于轴对称的特性，任意一点变形只发生在子午面上，因此可以作为 二 维问题处理。 6．等参数单元指的是：描述位移和描述坐标采用相同的形函数形式。等参数单元优点是：可以采用高阶次位移模式，能够模拟复杂几何边界，方便单元刚度矩阵和等效节点载荷的积分运算。 7．有限单元法首先求出的解是 节点位移 ，单元应力可由它求得，其计算公式为 {}{} [][]e D B σδ=。（用符号表示即可） 8．一个空间块体单元的节点有 3 个节点位移： u ，v ，w 9．变形体基本变量有位移应变应力 基本方程 平衡方程 物理方程 几何方程 10.实现有限元分析标准化和规范化的载体就是单元 《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0, 解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志 ,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n 《计算方法》2008试题与答案 一、填空题（每空2分，共20分） (1) 为了提高数值计算精度, 当正数x 充分大时, 应将)1ln(2--x x 改写为 _ln(x -______. (2) 3*x 的相对误差约是*x 的相对误差的_1/3____ 倍 (3).设?? ?? ? ?????---=283012251A ，则∞A =__13______.1A =___14_____ (4) 已知()p x 为二次多项式，满足(2)(2)3P f -=-=, (1)(1)1P f -=-=和 '(1)'(1)1P f -=-=，则()(2)(2)(2)(1)p x f a x b x x =-+++++，这里 a = -2 , b = 3 。 (5) 设32()4321f x x x x =+++，则差商[]3 ,2 ,1 ,0f =__4__[]0, 1, 2, 3, 4f =_0_. (6)n 个求积节点的求积公式的代数精确度最高为_21n -_____次. (7) 求解初值问题1)0(),(50'=+-=y x y y 时，若用改进欧拉方法的绝对稳定域中步长h 不 超过.0.04 二、（10分）用Newton 法求方程2ln =-x x 在区间) ,2(∞内的根, 取03x =， 要求 8110k k k x x x -+-<，计算过程中数值保留8位有效数字。 解 此方程在区间(2, )∞内只有一个根s ，而且在区间（2，4）内。设 ()ln 2f x x x =-- 则 ' 1()1f x x =- ， '' 21()f x x = Newton 法迭代公式为 1ln 2(1ln )11/1 k k k k k k k k x x x x x x x x +--+=- =--， （5分） 参考答案 第一章 1 *1x =1.7； * 2x =1.73； *3x =1.732 。 2． 3． （1） ≤++)(* 3*2*1x x x e r 0.00050； （注意：应该用相对误差的定义去求） （2） ≤)(*3*2*1x x x e r 0.50517； （3） ≤)/(*4*2x x e r 0.50002。 4．设6有n 位有效数字，由6≈2.4494……，知6的第一位有效数字1a =2。 令3)1()1(1* 102 1 102211021)(-----?≤??=?= n n r a x ε 可求得满足上述不等式的最小正整数n =4，即至少取四位有效数字，故满足精度要求可取6≈2.449。 5. 答：（1）*x (0>x )的相对误差约是* x 的相对误差的1/2倍； （2）n x )(* 的相对误差约是* x 的相对误差的n 倍。 6． 根据******************** sin 21)(cos 21sin 21)(sin 21sin 21)(sin 21)(c b a c e c b a c b a b e c a c b a a e c b S e r ++≤ =* *****) ()()(tgc c e b b e a a e ++ 注意当20* π < 7．设20= y ，41.1*0 =y ，δ=?≤--2* 00102 1y y 由 δ1* 001*111010--≤-=-y y y y ， δ2*111*221010--≤-=-y y y y M δ10*991*10101010--≤-=-y y y y 即当0y 有初始误差δ时，10y 的绝对误差的绝对值将减小10 10-倍。而110 10 <<-δ，故计算过程稳定。 8. 变形后的表达式为： （1）)1ln(2--x x =)1ln(2-+-x x （2）arctgx x arctg -+)1(=) 1(11 ++x x arctg （3） 1ln )1ln()1(ln 1 --++=? +N N N N dx x N N =ΛΛ+-+- +3 2413121)1ln(N N N N 1ln )11ln()1(-++ +=N N N N =1)1ln()1 1ln(-+++N N N （4）x x sin cos 1-=x x cos 1sin +=2x tg 计算方法上机报告 姓名: 学号: 班级: 目录 题目一----------------------------------------------------------------- 4 - 1.1题目内容-------------------------------------------------------- 4 - 1.2算法思想-------------------------------------------------------- 4 - 1.3Matlab 源程序----------------------------------------------------- 5 - 1.4计算结果及总结------------------------------------------------- 5 - 题目二----------------------------------------------------------------- 7 - 2.1题目内容-------------------------------------------------------- 7 - 2.2算法思想-------------------------------------------------------- 7 - 2.3 Matlab 源程序---------------------------------------------------- 8 - 2.4计算结果及总结------------------------------------------------- 9 - 题目三--------------------------------------------------------------- -11- 3.1题目内容----------------------------------------------------------- 11 - 3.2算法思想----------------------------------------------------------- 11 - 3.3Matlab 源程序--------------------------------------------------- -13 - 3.4计算结果及总结----------------------------------------------------- 14 - 题目四--------------------------------------------------------------- -15 - 4.1题目内容----------------------------------------------------------- 15 - 4.2算法思想----------------------------------------------------------- 15 - 4.3Matlab 源程序--------------------------------------------------- -15 - 4.4计算结果及总结----------------------------------------------------- 16 - 题目五--------------------------------------------------------------- -18 - -18 - 5.1题目内容 5.2算法思想----------------------------------------------------------- 18 - 5.3 Matlab 源程序--------------------------------------------------- -18 - 计算方法 一、填空题 1.假定x ≤1，用泰勒多项式?+??+++=! !212n x x x e n x ,计算e x 的值，若要求截断误差不超过0.005，则n=_5___ 2. 解 方 程 03432 3=-+x － x x 的牛顿迭代公式 )463/()343(121121311+--+--=------k k k k k k k x x x x x x x 3.一阶常微分方程初值问题 ?????= ='y x y y x f y 0 0)() ,(,其改进的欧拉方法格式为)],(),([21 1 1 y x y x y y i i i i i i f f h +++++= 4.解三对角线方程组的计算方法称为追赶法或回代法 5. 数值求解初值问题的四阶龙格——库塔公式的局部截断误差为o(h 5 ) 6.在ALGOL 中，简单算术表达式y x 3 + 的写法为x+y ↑3 7.循环语句分为离散型循环,步长型循环,当型循环. 8.函数)(x f 在[a,b]上的一次（线性）插值函数= )(x l )()(b f a b a x a f b a b x --+-- 9.在实际进行插值时插值时，将插值范围分为若干段，然后在每个分段上使用低阶插值————如线性插值和抛物插值，这就是所谓分段插值法 10、数值计算中，误差主要来源于模型误差、观测误差、截断误差和舍入误差。 11、电子计算机的结构大体上可分为输入设备 、 存储器、运算器、控制器、 输出设备 五个主要部分。 12、算式2 cos sin 2x x x +在ALGOL 中写为))2cos()(sin(2↑+↑x x x 。 13、ALGOL 算法语言的基本符号分为 字母 、 数字 、 逻辑值、 定义符四大 2020年奥鹏吉大网络教育《计算方法》大作业解答 (说明：前面是题目，后面几页是答案完整解答部分，注意的顺序。) 一、解线性方程 用矩阵的LU分解算法求解线性方程组 用矩阵的Doolittle分解算法求解线性方程组 用矩阵的Doolittle分解算法求解线性方程组 用高斯消去法求解线性方程组 用高斯消去法求解线性方程组 用主元素消元法求解线性方程组 用高斯消去法求解线性方程组 利用Doolittle分解法解方程组Ax=b，即解方程组 1、用矩阵的LU分解算法求解线性方程组 X1+2X2+3X3 = 0 2X1+2X2+8X3 = -4 -3X1-10X2-2X3 = -11 2、用矩阵的Doolittle分解算法求解线性方程组 X1+2X2+3X3 = 1 2X1– X2+9X3 = 0 -3X1+ 4X2+9X3 = 1 3、用矩阵的Doolittle分解算法求解线性方程组 2X1+X2+X3 = 4 6X1+4X2+5X3 =15 4X1+3X2+6X3 = 13 4、用高斯消去法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 5、用无回代过程消元法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 6、用主元素消元法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 7、用高斯消去法求解线性方程组 123123123234 4272266 x x x x x x x x x -+=++=-++= 8、利用Doolittle 分解法解方程组Ax=b ，即解方程组 12341231521917334319174262113x x x x -? ????? ???? ??-??????=? ? ????--?????? --???? ?? 计算方法上机报告 姓名： 学号： 班级： 目录 题目一------------------------------------------------------------------------------------------ - 4 - 1.1题目内容 ---------------------------------------------------------------------------- - 4 - 1.2算法思想 ---------------------------------------------------------------------------- - 4 - 1.3Matlab源程序----------------------------------------------------------------------- - 5 - 1.4计算结果及总结 ------------------------------------------------------------------- - 5 - 题目二------------------------------------------------------------------------------------------ - 7 - 2.1题目内容 ---------------------------------------------------------------------------- - 7 - 2.2算法思想 ---------------------------------------------------------------------------- - 7 - 2.3 Matlab源程序---------------------------------------------------------------------- - 8 - 2.4计算结果及总结 ------------------------------------------------------------------- - 9 - 题目三----------------------------------------------------------------------------------------- - 11 - 3.1题目内容 --------------------------------------------------------------------------- - 11 - 3.2算法思想 --------------------------------------------------------------------------- - 11 - 3.3Matlab源程序---------------------------------------------------------------------- - 13 - 3.4计算结果及总结 ------------------------------------------------------------------ - 14 - 题目四----------------------------------------------------------------------------------------- - 15 - 4.1题目内容 --------------------------------------------------------------------------- - 15 - 4.2算法思想 --------------------------------------------------------------------------- - 15 - 4.3Matlab源程序---------------------------------------------------------------------- - 15 - 4.4计算结果及总结 ------------------------------------------------------------------ - 16 - 题目五----------------------------------------------------------------------------------------- - 18 - 计算方法实验报告 班级： 学号： 姓名： 成绩： 1 舍入误差及稳定性 一、实验目的 （1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性 二、实验内容 1、用两种不同的顺序计算10000 21n n -=∑，分析其误差的变化 2、已知连分数() 1 01223//(.../)n n a f b b a b a a b =+ +++，利用下面的算法计算f ： 1 1 ,i n n i i i a d b d b d ++==+ (1,2,...,0 i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分 1 041 n n x y dx x =+? (0,1,...,1 n = 4、设2 2 11N N j S j == -∑ ，已知其精确值为1311221N N ?? -- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序 （3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数 三、实验步骤、程序设计、实验结果及分析 1、用两种不同的顺序计算10000 2 1n n -=∑，分析其误差的变化 （1）实验步骤： 分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算 #include 北京科技大学2004年《计算方法》试题及答案 一、判断题(下列各题，你认为正确的，请在括号内打“√”，错的打“×”,每题2分，共12分) 1、任何近似值的绝对误差总是大于其相对误差 （×） 2、3步Adams 隐式法比4步Adams 显式法的绝对稳定性要好。 （√） 3、在任何情况下，求解线性方程组时，Sidel 迭代法总是优于Jacobi 迭代法。 （×） 4、设],[)(b a C x f n ∈，若0)() (≡x f n ，],[b a x ∈，则0],,,[10=n x x x f ，其中 ],[b a x i ∈，n i ,,1,0 = （√） 5、给定n 个数据点，则至多构照1-n 次最小二乘多项式 （√） 6、数值求积公式的代数精确度越高，计算结果越可靠。 （×） 二、填空题(1、2、3小题每空1分，其他题每空2分,共20分) 1、设A 是一个108?的矩阵，B 是一个5010?的矩阵，C 是一个150?的矩阵，D 是 一个801?的矩阵,根据矩阵乘法结合率，ABCD F =可按如下公式计算 （1）[]D BC A F )(= （2）[])(CD B A F = 则公式（1）效率更高，其计算量为1240flops 。 2、设数据21,x x 的相对误差限分别为05.0和005.0，那么两数之商 2 1 x x 的相对误差限为 =)( 2 1 x x r ε0.055。 3、 设?? ? ???-=1123A ，则=1A 4，=∞A 5，=F A 15，=)(A ρ4，=∞)(A cond 4。 4、计算3 a 的割线法迭代公式为2 1 121 113133 1 )()(------++++=---=k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x 5、求解初值问题???=-='0)0() exp(2y x y 的改进后的 Euler 公式为 )]exp()[exp(2 2 121++-+-+=n n n n x x h y y 。 6、将正定矩阵???? ??????--=201032124A 作T LL 分解，则=L ?????? ???????? -8134 22 102 1 002 《计算方法B》上机 实验报告 学院：机械工程学院 班级: 姓名： 学号： 2015年12月22日 1 1.计算以下和式： S = ∑ 8n + 1 - 8n + 4 - 8n + 5 - 8n + 6 ? ，要求： 4 2 1 1 ∞ n =0 16 n ? ? ? ? (1)若保留 11 个有效数字，给出计算结果，并评价计算的算法； (2)若要保留 30 个有效数字，则又将如何进行计算。 实现思想： 以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和 以及被减数之和，最后将两者相减。另外，减数与被减数求和均为同号计算， 按照绝对值递增顺序相加可减小舍入误差。此题中对有效数字有要求，因而计 算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。 源程序： m=input('输入有效数字个数m='); s0=1;s1=0;s2=0;n=0; %判断迭代次数 while s0>=0.5*10^-(m-1) s0=4/(16^n*(8*n+1))-2/(16^n*(8*n+4))-1/(16^n*(8*n+5))- 1/(16^n*(8*n+6)); n=n+1; end %分别求解各项并求和 for k=n-1:-1:0 a1=4/(16^k*(8*k+1)); a2=2/(16^k*(8*k+4)); a3=1/(16^k*(8*k+5)); a4=1/(16^k*(8*k+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=vpa(s1-s2,m) 实验结果：11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。 图1.11位有效数字计算结果图2.30为有效数字计算结果 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 北京市科技计划项目课题经费预算编制说明 一、项目课题建议单位在申报项目课题时必须按照北京市科委的规定编制项目课题经费预算。项目课题经费预算的编制应严格遵守目标相关、政策相符、实事求是的原则。 二、课题经费来源包括项目课题建议单位从各个渠道筹集的资金： 市科技经费：指北京市科委拨付的财政经费，包括科三费和科学事业费。 国家有关部委拨款：指除北京市科委以外的其它各级政府（如国家科技部、北京市计委、区科委等）为实施本项目课题拨付的财政经费。 项目依托单位匹配经费：指具有直接行政隶属的主管单位为实施本项目课题拨付的经费。 课题承担单位自筹经费：企业通过其它渠道筹措到的资金，如股东增资、历年的经营利润等。 银行贷款：指为实施本项目课题，项目课题建议单位从经中国人民银行批准的可以经营信贷业务的金融机构处获得的贷款经费。 其它：指未列入以上各项的其它经费来源。 三、课题经费支出即项目课题研究过程中发生的全部费用支出预算： 人员费用：指直接参加项目课题研究人员的工资性费用，包括专职人员费用及外聘人员费用。列入的人员要与项目课题任务书中参加的人员一致，其中：项目课题组人员所在单位有事业费拨款的，由所在单位按照国家规定的标准从事业费中及时足额支付给项目课题组成员，并按规定在项目课题预算的相关科目中列示，不得在国家资助的项目课题经费中重复列支。国家另有规定的，按照有关规定执行。 试验外协费：指研究、开发项目课题所发生的带料外加工或因本单位不具备条件面委托外单位进行试验、加工、测试、计算等发生的费用。发生试验外协费时，必须与协作单位签订合同书。 合作交流费：指项目课题研究过程中需与国内外机构开展合作研究所发生的费用。发生合作费时，必须与合作机构签订相关的合同书。 数值计算第一次大作业 实验目的 以Hilbert 矩阵为例，研究处理病态问题可能遇到的困难。 内容 Hilbert 矩阵的定义是 ,() 11/21/31/1/21/31/41/(1)1/31/41/51/(2)1/1/(1)1/(2)1/(21)n i j H h n n n n n n n =????+????=+??????++-?? 它是一个对称正定矩阵，而且()n cond H 随着n 的增加迅速增加，其逆矩阵1,()n i j H α-=， 这里 ,2(1)(1)!(1)!(1)[(1)!(1)!]()!()! i j i j n i n j i j i j n i n j α+-+-+-=+----- 1) 画出ln(())~n cond H n 之间的曲线（可以用任何的一种范数）。你能猜出 ln(())~n cond H n 之间有何种关系吗？提出你的猜想并想法验证。 用行范数 for n=1:50 for i=1:n for j=1:n A(i,j)=1/(i+j-1); B(i,j)=factorial(n+i-1)*factorial(n+j-1)/((i+j-1)*(factorial(i-1)*facto rial(j-1))^2*factorial(n-i)*factorial(n-j)); end end result1=0; for j=1:n result1=result1+A(1,j); end result1=log(result1); result2=0; for i=1:n for j=1:n result2=B(i,j)+result2; end result(i)=log(result2); end m=max(result); 第一次&第二次上机作业 上机作业： 1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5 给出执行结果，并简要分析一下产生现象的原因。 解：执行结果如下： 在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。 2.（课本181页第一题） 解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码 从以上代码显示的结果可以看出，I 20的近似值为0.7465 （2）I I =∫I I 5+I 10dx,可得∫I I 610dx ≤∫I I 5+I 10dx ≤∫I I 510dx,得 16(I +1)≤I I ≤15(I +1),则有1126≤I 20≤1105, 取I 20=1 105 ,以此逆序估算I 0。代码段及结果如下图所示 （3）从I20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设I I表示I I的近似值，I I-I I=(?5)I(I0?I0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。 2.（课本181页第二题） （1）上机代码如图所示 求得近似根为0.09058 （2）上机代码如图所示 得近似根为0.09064； （3）牛顿法上机代码如下 计算所得近似解为0.09091 第三次上机作业上机作业181页第四题 线性方程组为 [1.13483.8326 0.53011.7875 1.16513.4017 2.53301.5435 3.4129 4.9317 1.23714.9998 8.76431.3142 10.67210.0147 ][ I1 I2 I3 I4 ]=[ 9.5342 6.3941 18.4231 16.9237 ] （1）顺序消元法 A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435; 3.4129, 4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147]; b=[9.5342;6.3941;18.4231;16.9237]; 上机代码（函数部分）如下 function [b] = gaus( A,b )%用b返回方程组的解 B=[A,b]; n=length(b); RA=rank(A); RB=rank(B);北科大研究生计算方法作业
数值计算方法大作业
计算方法上机题答案
北京科技大学有限元试题及答案
《数值计算方法》上机实验报告
北京科技大学计算方法试题
西工大计算方法作业答案
西安交通大学计算方法B大作业资料
计算方法试题库讲解
2020年奥鹏吉大网络教育《计算方法》大作业解答
西安交通大学计算方法B大作业
计算方法上机实习题大作业(实验报告).
北京科技大学2004年《计算方法》试题及答案
西安交大计算方法b大作业课件
(完整版)数值计算方法上机实习题答案
北京科技大学 经费预算说明
计算方法大作业第一次
计算方法上机作业集合
- 计算方法大作业作业((北京科技大学研究生结课考试)
- 计算方法作业
- 数值计算方法大作业
- 计算流体力学课程大作业.doc
- 工程计算方法及软件应用--本科生考查大作业
- 西安交大计算方法B2017大作业
- 计算方法大作业——三次样条插值
- 解线性方程组的迭代法_计算方法大作业
- 西安交通大学计算方法B大作业
- 计算方法大作业----------雅可比
- 数值计算方法大作业
- 计算方法大作业——龙贝格积分法的matlab程序
- 数值计算方法大作业
- 数值计算方法大作业
- 计算方法大作业非线性方程求根的新方法
- 计算方法大作业
- 浙江大学计算方法大作业第十二题
- 《应用计算方法教程》大作业
- 2020年奥鹏吉大网络教育《计算方法》大作业解答
- 2018-2019学年第二学期期末考试《计算方法》大作业参考答案