计算方法作业第六章
1．考虑两个线性方程组，其系数矩阵如下
1211
11...23211111...1212341,121
1111...3452..............................121111...
12
21n n A A n n n n n ?
????
?-??????--??
+??
????==--??????+????
??-??
?
?????++-??
问题的真解均取为[1,1,1,1,...1]T x =，线性方程组的右端项用这个真解计算出来。相应的问题分别称为问题I 和问题II 。请进行如下数值实验：
（1） 对问题I 分别用Gauss 消元法，Cholesky 方法，修改的LDLT 算法，追赶法四
种方法求解，其中n=100；
（2） 对问题II 分别用Gauss 消去法，列主元Gauss 消去法，不做行交换的列主元
Gauss 消去法求解，其中n=6；
（3） 不断增加问题II 的矩阵阶数n=6,8,10,…,20，重复（2）的工作，看看会有什么
问题发生？解释其原因。
（1） Gauss ： 计算程序： n=100; A=2*eye(n); for i=1:n-1 A(i+1,i)=-1; A(i,i+1)=-1; end b=0; b(1)=1; b(100)=1;
[x,XA]=GaussJordanXQ(A,b); Gauss 消元法源程序：
%用Gauss 消元法解线性方程组 function [x,XA]=GaussJordanXQ(A,b) N = size(A); n = N(1); for i=1:(n-1)
for j=(i+1):n
if(A(i,i)==0)
disp('对角元素为0！');
%防止对角元素为0
return;
end
l = A(j,i);
m = A(i,i);
A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;
%消元方程
b(j)=b(j)-l*b(i)/m;
end
end
x=SolveUpTriangle(A,b);
%通用的求上三角系数矩阵线性方程组的函数XA = A;
%消元后的系数矩阵
（SolveUpTriangle.m）解上三角方程组源程序：%解上三角方程组
function x=SolveUpTriangle(A,b)
N = size(A);
n = N(1);
x(n)=b(n)/A(n,n);
for k=n-1:1
s=0;
for i=k+1:n
s=s+A(k,i)*x(i);
end
x(k)=(b(k)-s)/A(k,k);
end
结果：
x1=[0,0,0,…..0,0,1]T
x2=[0,0,0,…..0,0,0.3820]T
x3=[0,0,0,…..0,0,0.9900]T
x4=[1,1,1,…..1,1,1]T
Cholesky：
计算程序：
n=100;
A=2*eye(n);
for i=1:n-1
A(i+1,i)=-1;
A(i,i+1)=-1;
end
b=0;
b(1)=1;
b(100)=1;
x2=Cholesky(A,b); (Cholesky.m):Cholesky法源程序%用Cholesky方法解线性方程组function x=Cholesky(A,b);
N=size(A);
n=N(1);
L(1,1)=sqrt(A(1,1));
%L为下三角矩阵
for j=2:n
L(j,1)=A(1,j)./L(1,1)
end
for k=2:n
s=0;
for i=1:k-1
s=s+L(k,i);
end
L(k,k)=sqrt(A(k,k)-s);
for j=(k+1):n
sum=0;
for i=1:k-1
sum=sum+L(j,i)*L(k,i);
end
L(j,k)=(A(j,k)-sum)./L(k,k); end
end
%定义L的转置
for j=1:n
for k=1:n
LT(j,k)=L(k,j);
end
end
y=SolveDownTriangle(L,b);
x=SolveUpTriangle(LT,y);
解上三角方程组源程序：
%解上三角方程组
function x=SolveUpTriangle(A,b)
N = size(A);
n = N(1);
x(n)=b(n)/A(n,n);
for k=n-1:1
s=0;
for i=k+1:n
s=s+A(k,i)*x(i);
end
x(k)=(b(k)-s)/A(k,k);
End
解下三角方程组源程序：function x=SolveDownTriangle(A,b) N = size(A);
n = N(1);
x(1)=b(1)/A(1,1);
for k=2:n
s=0;
for i=1:k-1
s=s+A(k,i)*x(i);
end
x(k)=(b(k)-s)/A(k,k);
end
结果：
x1=[0,0,0,…..0,0,0.9894]T;
x2=[0,0,0,…..0,0,0.9899]T;
LTDT：
计算程序：
n=100;
A=2*eye(n);
for i=1:n-1
A(i+1,i)=-1;
A(i,i+1)=-1;
end
b=0;
b(1)=1;
b(100)=1;
x3=GJCholesky(A,b);
(GJCholesky.m):Cholesky法源程序
%用改进Cholesky法求解线性方程组
%记LD为U，U为下三角矩阵，LT为单位上三角矩阵function x=GJCholesky(A,b);
N=size(A);
n=N(1);
for k=1:n
U(k,1)=A(k,1);
end
for i=2:n
L(k,1)=A(1,k)./U(1,1);
end
for k=2:n
for j=k:n
s=0;
for i=1:k-1
s=s+U(j,i).*L(j,i);
end
U(j,k)=A(j,k)-s;
end
for j=k+1:n
sum=0;
for i=1:k-1
sum=sum+U(k,i).*L(j,i);
end
L(j,k)=(A(k,j)-sum)./U(k,k);
end
end
for k=1:n
D(k,k)=U(k,k);
L(k,k)=1;
LT(k,k)=1;
for j=k+1:n
LT(k,j)=L(j,k);
end
end
y=SolveDownTriangle(L,b);
for i=1:n
z(i)=y(i)./D(i,i);
end
x=SolveUpTriangle(LT,z);
解上三角方程组源程序：
%解上三角方程组
function x=SolveUpTriangle(A,b)
N = size(A);
n = N(1);
x(n)=b(n)/A(n,n);
for k=n-1:1
s=0;
for i=k+1:n
s=s+A(k,i)*x(i);
end
x(k)=(b(k)-s)/A(k,k);
end
解下三角方程组源程序：function x=SolveDownTriangle(A,b) N = size(A);
n = N(1);
x(1)=b(1)/A(1,1);
for k=2:n
s=0;
for i=1:k-1
s=s+A(k,i)*x(i);
end
x(k)=(b(k)-s)/A(k,k);
end
结果：
n=8
x1=[0,0,0,…..0,0,0.9891]T;x2=[0,0,0,…..0,0,0.9889]T; n=10
x1=[0,0,0,…..0,0,0.9889]T;x2=[0,0,0,…..0,0,0.9885]T; n=12
x1=[0,0,0,…..0,0,0.9886]T;x2=[0,0,0,…..0,0,0.9885]T; n=14
x1=[0,0,0,…..0,0,0.9884]T;x2=[0,0,0,…..0,0,0.9884]T; n=16
x1=[0,0,0,…..0,0,0.9881]T;x2=[0,0,0,…..0,0,0.9882]T; n=18
x1=[0,0,0,…..0,0,0.9878]T;x2=[0,0,0,…..0,0,0.9879]T; n=20
x1=[0,0,0,…..0,0,0.9875]T;x2=[0,0,0,…..0,0,0.9875]T;
追赶法：
计算程序：
n=100;
A=2*eye(n);
for i=1:n-1
A(i+1,i)=-1;
A(i,i+1)=-1;
end
b=0;
b(1)=1;
b(100)=1;
x4=zhuigan(A,b);
追赶法源程序
%用追赶法求解线性方程组
function M=zhuigan(A,g)
n=length(A);
L=eye(n);
U=zeros(n);
for i=1:n-1
U(i,i+1)=A(i,i+1);
end
U(1,1)=A(1,1);
for i=2:n
L(i,i-1)=A(i,i-1)/U(i-1,i-1);
U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);
end
Y(1)=g(1);
for i=2:n
Y(i)=g(i)-L(i,i-1)*Y(i-1);
end
M(n)=Y(n)/U(n,n);
for i=n-1:-1:1
M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i); end
（2）
Gauss消去：
计算程序：
n=6;
for i=1:n
s=0;
for j=1:n
A(i,j)=1/(i+j-1);
s=s+A(i,j);
end
b(i)=s;
end
[x1,XA]=GaussJordanXQ(A,b); Gauss消元法源程序
%用Gauss消元法解线性方程组function [x,XA]=GaussJordanXQ(A,b) N = size(A);
n = N(1);
for i=1:(n-1)
for j=(i+1):n
if(A(i,i)==0)
disp('对角元素为0！');
%防止对角元素为0
return;
end
l = A(j,i);
m = A(i,i);
A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;
%消元方程
b(j)=b(j)-l*b(i)/m;
end
end
x=SolveUpTriangle(A,b);
%通用的求上三角系数矩阵线性方程组的函数XA = A;
%消元后的系数矩阵
解上三角方程组源程序：
%解上三角方程组
function x=SolveUpTriangle(A,b)
N = size(A);
n = N(1);
x(n)=b(n)/A(n,n);
for k=n-1:1
s=0;
for i=k+1:n
s=s+A(k,i)*x(i);
end
x(k)=(b(k)-s)/A(k,k);
end
列主元Gauss消去：
计算程序：
n=6;
for i=1:n
s=0;
for j=1:n
A(i,j)=1/(i+j-1);
s=s+A(i,j);
end
b(i)=s;
end
[x2,XA]=GaussLineMainXQ(A,b);
列主元Gauss消去法源程序
%列主元Gauss消去法解线性方程组function x=GaussLinemainXQ(A,b); N=size(A);
n=N(1);
index=0;
for i=1:n-1
main=max(abs(A(1:n,i)));
%选取列主元
for k=i:n
if (abs(A(k,i)==main)
index=k;
break;
end
%保存主元所在行的指标
end
temp=A(i,1:n);
A(i,1:n)=A(index,1:n);
A(index,1:n)=temp;
btemp=b(index);
b(index)=b(i);
b(i)=btemp;
%交换主元所在行
for j=(i+1):n
if (A(i,i)==0)
disp('对角线元素为0’);
return;
end
l=A(j.i);
m=A(i,i);
A(j,1:n)=A(j,1:n)-A(i,1:n).*l/m; b(j)=b(j)-b(i).*l/m;
end
end
x=SolveUpTriangle(A,b);
解上三角方程组源程序：
%解上三角方程组
function x=SolveUpTriangle(A,b) N = size(A); n = N(1); x(n)=b(n)/A(n,n); for k=n-1:1 s=0; for i=k+1:n s=s+A(k,i)*x(i); end
x(k)=(b(k)-s)/A(k,k); end
2．考虑矩阵A I B =+，其中(),ij ij B b b i j ==-。利用Gauss-Jordan 方法求解逆矩阵，并验证逆矩阵与A 的乘积是否为单位矩阵 Gauss-Jordan 方法程序：
function x=GaussJordan(A) [n,n]=size(A); while det(A)==0 break end
s=1;P=zeros(1,n); while s<=n
max=abs(A(s,s)); big=0; for m=s:n
if max if big==1 P(s)=q; D=A(s,:); A(s,:)=A(q,:); A(q,:)=D; end for i=1:n if i~=s for j=1:n if j~=s A(i,j)=A(i,j)-A(s,j)*A(i,s)/A(s,s); else continue end end else continue end A(i,s)=-A(i,s)/A(s,s); end for k=1:n if k A(s,k)=A(s,k)/A(s,s); end if k==s A(s,k)=1/A(s,s); end if k>s A(s,k)=A(s,k)*A(s,s); end end s=s+1; end for i=n:(-1):1 if P(i)~=0 p=P(i); d=A(:,i); A(:,i)=A(:,p); A(:,p)=d; else continue end end A 调用程序为： format rat n=6; a=zeros([n,n]); for i=1:n for j=1:n a(i,j)=i-j; End end b=eye([n,n])+a B=GaussJordan(b) B*b 结果为： B = 51/106 -39/106 -23/106 -7/106 9/106 25/106 -41/106 75/106 -21/106 -11/106 -1/106 9/106 -27/106 -23/106 87/106 -15/106 -11/106 -7/106 -13/106 -15/106 -17/106 87/106 -21/106 -23/106 1/106 -7/106 -15/106 -23/106 75/106 -39/106 15/106 1/106 -13/106 -27/106 -41/106 51/106 B*b ans = 1.0000 0.0000 0 -0.0000 -0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 1.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 是单位矩阵 第六章运算方法 1 下列数中最小的数为——。 A．(101001)2 B (52)8 C （2B）16 2．下列数中最大的数为。 A．(10010101)2 B．(227)d C．(96)16 3．设寄存器位数为8位，机器数采用补码形式(含1位符号位)，对应于十进制数（-27），寄存器内容为一——。 A．27H B．9BH C．E5K 4．对真值0表示形式唯一的机器数是——o A．原码B．补码和移码 C 反码 D 以上都不对 5． 6 在整数定点机中，下述正确的说法是 A．原码和反码不能表示—1，补码可以表示—1 B．三种机器数均可表示—1 c．三种机器数均可表示—1，且三种机器数的表示范围相同 7在小数定点机中，下述说法正确的是——。 A．只有补码能表示—1 B．只有原码不能表示—1 c．三种机器数均不能表示—1 8．某机字长8位．采用形式(其中1位为符号位)则机器数所能表示的范围 A．一127—127 D．一128，十128 C 一128一十127 9、用n+1位字长表示定点数(其中1位为符号位)，它所能表示的整数范围是 能表示的小数范围是。 A、阶码取4位(台阶符1位)，尾数取12位(合数符1位) B．阶码取5位(台阶符1位)，尾数取11位(合数符1位) c．阶码取8位(含阶符1位)，尾数取8位(合数符1位) 70在下述有关不恢复余数法何时需恢复余数的说法中，——是正确的A最后一次余数为正时，要恢复 B．最后一次余数为负时，要恢复 C．最后一次余数为。时，要恢复 D．任何时候都不恢复余数 71．在定点机中执行算术运其时会产生溢出，其原因是——。 A．主存容量不够B．运算结果无法表示 c．操作数地址过大D．以上都对 72．在浮点机中，下列说法是正确的。 A．尾数的第一数位为1时，即为规格化形式 B、尾数的第一数值与数符不同时，即为规格化形2 目录 第一章非线性方程求根 (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） 计算方法上机报告 姓名： 学号： 班级： 上课班级： 说明： 本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。 1. 对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=0681581482184161n n n n S n ，要求： ① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算； （1） 算法思想 1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114 16818485861681 n n n a 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 10m t -≤ end; 3. for ,1,2,,0n i i i =--??? ;s s t =+ （3）Matlab源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0; for n=0:50 t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); if t<=10^(-m) %判断通项与精度的关系break; end end; fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值 for i=n-1:-1:0 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s=s+t; %求和运算 end s=vpa(s,m) %控制s的精度 （4）结果与分析 当保留11位有效数字时，需要将n值加到n=7， s =3.1415926536； 当保留30位有效数字时，需要将n值加到n=22, s =3.14159265358979323846264338328。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 参考答案 第一章 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 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 -? ????? ???? ??-??????=? ? ????--?????? --???? ?? 94 第六章 计算方法简介 §1 数值逼近 1.1 插值 许多实际问题都要用函数)(x f y =来表示某种内在规律的数量关系，其中相当一部分函数虽然可能在某个区间上具有很好的性质（连续、光滑等），但没有函数的表达式信息，我们只能通过实验或者观测得到函数在一些点i x 上的函数值 )(i i x f y =),2,1,0(n i =，这是一张函数表.有些函数虽然有解析式，但由于计算 复杂，使用不方便，我们通常也造一个函数表，例如三角函数表、平方根表等. 为了研究函数的性质，往往还需要求出不在函数表上的函数值，因此我们希望根据给定的函数表构造一个既能反映函数)(x f y =的性质、又便于计算的简单函数 )(x P ，用)(x P 来近似)(x f .这就是插值所要研究的问题. )(x P 称为)(x f 的插值函数.常用的插值函数是代数多项式或分段代数多项式. 1.1 Lagrange 插值 1.1.1 方法介绍 Lagrange 插值方法即，给定n 个插值节点以及对应的函数值信息， )(i i x f y =),2,1,0(n i =，利用n 次Lagrange 插值多项式公式，则对插值区间内 任意x 的函数值y 可通过下式近似求得： )()(1 1 ∏ ∑≠==--=n k j j j k j n k k x x x x y x y . 其中 ∏≠=--n k j j j k j x x x x 1称为插值基函数.可见，在Lagrange 插值中，对应1+n 个节点的 插值基函数一共有1+n 个，每个基函数是一个n 次多项式. 1.1.2 MATLAB 实现 Lagrange.m 计算方法上机报告 姓名: 学号: 班级: 目录 题目一----------------------------------------------------------------- 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 舍入误差及稳定性 一、实验目的 （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 《计算方法》上机指导书 实验1 MATLAB 基本命令 1．掌握MATLAB 的程序设计 实验内容：对以下问题，编写M 文件。 (1) 生成一个5×5矩阵，编程求其最大值及其所处的位置。 (2) 编程求∑=20 1!n n 。 (3) 一球从100米高度自由落下，每次落地后反跳回原高度的一半，再落下。求它在 第10次落地时，共经过多少米？第10次反弹有多高？ 2．掌握MATLAB 的绘图命令 实验内容：对于自变量x 的取值属于[0，3π]，在同一图形窗口画出如下图形。 （1）1sin()cos()y x x =?； （2）21 2sin()cos()3 y x x =-； 实验2 插值方法与数值积分 1. 研究人口数据的插值与预测 实验内容：下表给出了从1940年到1990年的美国人口，用插值方法推测1930年、1965年、2010年人口的近似值。 美国人口数据 1930年美国的人口大约是123,203千人，你认为你得到的1965年和2010年的人口数字精确度如何？ 2．最小二乘法拟合经验公式 实验内容：某类疾病发病率为y ‰和年龄段x (每五年为一段，例如0~5岁为第一段，6~10岁为第二段……)之间有形如bx ae y =的经验关系，观测得到的数据表如下 （1）用最小二乘法确定模型bx ae y =中的参数a 和b 。 （2）利用MATLAB 画出离散数据及拟合函数bx ae y =图形。 3. 复化求积公式 实验内容：对于定积分? +=1 02 4dx x x I 。 （1）分别取利用复化梯形公式计算，并与真值比较。再画出计算误差与n 之间的曲线。 （2）取[0，1]上的9个点，分别用复化梯形公式和复化辛普森公式计算，并比较精度。 计算方法大作业 题目：非线性方程求根的新方法 班级：xxx 学号：xxx 姓名：xxx 非线性方程求根的新方法 一、问题引入 在计算和实际问题中经常遇到如下非线性问题的求解： F(x)=0 (1) 我们经常采用的方法是经典迭代法： 经典迭代方法 不动点迭代方法是一种应用广泛的方法，其加速方法较多，如Stiffensen加速方法的局部收敛阶(以下简称为收敛阶)为2阶；牛顿迭代方法的收敛阶亦为2阶，且与其相联系的一些方法如简化牛顿法、牛顿下山法、弦截法的收敛阶阶数介于1和2之间；而密勒法的收敛阶与牛顿法接近，但计算量较大且涉及零点的选择问题，同时收敛阶也不够理想。 因此本文介绍一种新的迭代方法 从代数角度看，牛顿法和密勒法分别是将f(x)在xk附近近似为一线性函数和二次抛物插值函数，一种很自然的想法就是能否利用Taylor展开，将f(x)在xk附近近似为其他的二次函数?答案是肯定的.其中的一种方法是将f(x)在Xk处展开3项，此时收敛阶应高于牛顿法，这正是本文的出发点． 二、算法推导 设函数f(x)在xk附近具有二阶连续导数，则可将f(x)在xk处进行二阶Taylor展开，方程(1) 可近似为如下二次方程： f(xk)+f’(xk)(x-xk)+2^(-1)f’’(xk)(x-xk)^2=0，(2) 即 2^(-1)f’’(xk)x^2+(f’(xk)-xkf’’(xk))x+2^(-1)f’’(xk)xk^2-xkf’(xk)+f(xk)=0(3) 利用求根公式可得 X=xk-(f’’(xk))^(-1)(f’(xk))-sqrt((f’(xk)^2±2f’’(xk)f(xk)))(4) 其中±符号的选取视具体问题而定，从而可构造迭代公式 X k+1=xk-(f’’(xk))^(-1)(f’(xk))-sqrt((f’(xk)^2±2f’’(xk)f(xk)))(5) 确定了根号前正负号的迭代公式(5)，可称为基于牛顿法和Taylor展开的方法，简记为BNT 方法． 为描述方便起见，以下将f(xk)，f’(xk)，f’’(xk)分别记为f,f’,f’’．首先，二次方程(3)对应于一条抛物曲线，其开口方向由f’’(xk)，x∈U(xk)的符号确定，其中U(xk)为xk的某邻域，其顶点为 P(xk-(f’’)^(-1)f’，fk-(2f’’)^(-1)(f’)^2)．为使(5)式唯一确定x k+1，须讨论根式前正负号的取舍问题．下面从该方法的几何意义分析(5)式中正负号的取舍． 1)当f(xk)=o时，z。即为所求的根． 2)当f(xk)>O时，根据y=f(x)的如下4种不同情形(见图1)确定(5)式中根号前的符号． (a)当f’’(xk) 工程计算方法与软件应用 本科生大作业 考核方式：考查（成绩按各软件的课外作业成绩综合给出）。 各软件讲完后1~2星期内上交作业。 一、CAD/CAE软件作业（每个学生完成下列任意一题） 题目一： 一端固定支撑，一端集中力的梁，横截面为10x10cm，长为150cm，受集中载荷作用，P=50N。弹性模量E=70GPa，泊松比r=0.2。用ABAQUS 软件建模并计算最大应力和最大位移的位置和大小。 （1）二维；（2）三维 图1梁受力简图 题目二： 图中所示为一个连接件，一端焊接到设备母体上，一端在圆柱销子作用下的圆孔，圆孔下半周受到30 kN的均布载荷作用，用ABAQUS 软件建模并计算最大应力和最大位移的位置和大小。 图2 连接件受力简图 题目三： 如图3所示为一薄壁圆筒，在圆筒中心受集中力F作用，对此进行受力分析，并给出应力、位移云图，并求A、B两点位移。 圆筒几何参数：长度L=0.2m；半径R=0.05m壁厚t=2.5mm。 材料参数：弹性模量E=120Gpa；泊松比0.3 载荷：F=1.5kN。 图3薄壁管受力简图 题目四： 如图4所示为一燃气输送管道截面及受力见图，试分析管道在内部压力作用下的应力场。 几何参数：外径0.6m，内径0.4m，壁厚0.2m 材料参数：弹性模量E=120Gpa；泊松比0.26 载荷P=1Mpa。 图4燃气管受力简图 题目五： 如图5为一三角桁架受力简图，途中各杆件通过铰链链接，杆件材料及几何参数见表1和表2所示，桁架受集中力F1=5kN、F2=2.5kN 作用，求桁架各点位移及反作用力。 图5 三角桁架受力简图 表1 杆件材料参数 表2 杆件几何参数 第六章 数值计算方法举例 6-1二分法(Bisection) 二分法是最简单的解法，该算法只有简单的几个步骤。 （1）先猜两个值a 、b ，使得f(a)*f(b)小于0，也就是f(a)、f(b)必须异号。这样才能保证在a 与b 间存在一个c 值，使得f(c)=0。 （2）令c=(a+b)/2,如果f(c)=0就找到了一个解，工作完成。 （3）f(c)不为时，如果f(a)、f(b)异号，则以a 、c 为新的两个猜测值来重复步骤2；如果f(b)、f(c)异号，则以b 、c 为新的猜测值来重复步骤2。 编程举例：用二分法求解方程sin x x =及方程2200x x --=的根。 1. program main 2. implicit none 3. real :: a, b !两个猜测值 4. real :: ans !算出的值 5. real, external :: bisect, f1, f2 6. do while(.true.) 7. write(*,*)＂输入两个猜测值＂ 8. read(*,*)a, b 9. !f(a)*f(b)<0的猜测值才是有效的猜测 10. if(f1(a)*f1(b)<0)exit 11. write(*,*)＂不正确的猜测＂ 12. end do 13. !调用二分法求根的函数 14. ans=bisect(a, b, f1) 15. ！显示结果 16. write(*,＂(＇x=＇,F6.3)＂) ans 17. 18. do while(.true.) 19. write(*,*)＂输入两个猜测值＂ 20. read(*,*)a, b 21. !f(a)*f(b)<0的猜测值才是有效的猜测 22. if(f2(a)*f2(b)<0)exit 23. write(*,*)＂不正确的猜测＂ 24. end do 25. !调用二分法求根的函数 26. ans=bisect(a, b, f2) 27. ！显示结果 28. write(*,＂(＇x=＇,F6.3)＂) ans 29. stop 30. end 31. real function bisect(a, b, func) 32. implicit none 计算方法大作业 机械电子工程系 老师：廖福成 注：本文本只有程序题，证明题全部在手写已交到理化楼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 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 计算方法A 上机大作业 1. 共轭梯度法求解线性方程组 算法原理：由定理3.4.1可知系数矩阵A 是对称正定矩阵的线性方程组Ax=b 的解与求解二次函数1()2 T T f x x Ax b x =-极小点具有等价性，所以可以利用共轭梯度法求解1()2 T T f x x Ax b x = -的极小点来达到求解Ax=b 的目的。 共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式： (1)()()k k k k x x d α+=+ 产生的迭代序列(1)(2)(3)x x x ，，，... 在无舍入误差假定下，最多经过n 次迭代，就可求得()f x 的最小值，也就是方程Ax=b 的解。 首先导出最佳步长k α的计算式。 假设迭代点()k x 和搜索方向()k d 已经给定，便可以通过()()()() k k f x d φαα=+的极小化 ()()min ()()k k f x d φαα=+ 来求得，根据多元复合函数的求导法则得: ()()()'()()k k T k f x d d φαα=?+ 令'()0φα=，得到: ()() ()()k T k k k T k r d d Ad α=，其中()()k k r b Ax =- 然后确定搜索方向()k d 。给定初始向量(0)x 后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0) (0)(0)(0)()d r f x b Ax ==-?=-。令 (1)(0)00x x d α=+ 其中(0)(0)0(0)(0) T T r d d Ad α=。第二次迭代时，从(1) x 出发的搜索方向不再取(1)r ，而是选取(1) (1)(0)0d r d β=+，使得(1)d 与(0)d 是关于矩阵A 的共轭向量，由此可 求得参数0β: 第一次&第二次上机作业 上机作业： 1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5 给出执行结果，并简要分析一下产生现象的原因。 解：执行结果如下： 在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。 2.（课本181页第一题） 解：（1）n=0时，积分得=ln6-ln5，编写如下图代码 从以上代码显示的结果可以看出，的近似值为0.7465 （2）=dx,可得dx dx dx,得 ,则有, 取,以此逆序估算。代码段及结果如下图所示 结果是从逆序输出至，所以得到的近似值为0.0227。（3）从估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设表示的近似值，-=(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。 2.（课本181页第二题）（1）上机代码如图所示 求得近似根为0.09058 （2）上机代码如图所示 得近似根为0.09064； （3）牛顿法上机代码如下 计算所得近似解为0.09091 第三次上机作业 上机作业181页第四题 线性方程组为 = （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返回方程组的解 题目利用数值计算方法求取基尼系数 姓名与学号 指导教师 年级与专业 所在学院 一、问题综述： 基尼系数（Gini coefficient），是20世纪初意大利学者科拉多·吉尼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。是比例数值，在0和1之间。基尼指数（Gini index）是指基尼系数乘100倍作百分比表示。在民众收入中，如基尼系数最大为“1”，最小等于“0”。前者表示居民之间的收入分配绝对不平均（即所有收入都集中在一个人手里，其余的国民没有收入），而后者则表示居民之间的收入分配绝对平均，即人与人之间收入绝对平等，但这两种情况只出现在理论上；因此，基尼系数的实际数值只能介于0～1之间，基尼系数越小收入分配越平均，基尼系数越大收入分配越不平均。 设右图中的 实际收入分配曲线 （红线）和收入分 配绝对平等线（绿 线）之间的面积为 A，和收入分配绝 对不平等线（蓝 线）之间的面积为 B，则表示收入与 人口之间的比例的基尼系数为 A A+B 。 如果A为零，即基尼系数为0，表示收入分配完全平等（红线和绿线重叠）；如果B为零，则系数为1，收入分配绝对不平等（红线和蓝线重叠）。该系数可在0和1之间取任何值。实际上，一般国家的收入分配，既不是完全平等，也不是完全不平等，而是在两者之间，劳伦茨曲线为一条凸向横轴的曲线。收入分配越趋向平等，劳伦茨曲线的弧度越小（斜度越倾向45度），基尼系数也越小；反之，收入分配越趋向不平等，劳伦茨曲线的弧度越大，那么基尼系数也越大。 基尼系数的调节需要国家通过财政政策进行国民收入的二次分配，例如对民众的财政公共服务支出和税收等，从而让收入均等化，令基尼系数缩小。 基尼系数由于给出了反映居民之间贫富差异程度的数量界线，可以较客观、直观地反映和监测居民之间的贫富差距，预报、预警和防止居民之间出现贫富两极分化。因此得到世界各国的广泛认同和普遍采用。 联合国有关组织规定： ●若低于0.2表示收入平均； ●0.2-0.3表示相对平均； ●0.3-0.4表示相对合理； ●0.4-0.5表示收入差距大； ●0.6以上表示收入差距悬殊。 2013年1月18日，中国国家统计局一次性公布了自2003年以来十年的全国基尼系数。大陆统计局局长马建堂称，按照国际新的统计口径，大陆居民收入的基尼系数，2003年是0.479，2004年是0.473，2005年为0.485，2006年为0.487，2007年为0.484，2008年为0.491，2009年为0.490，2010年为 0.481，2011年为0.477，到2012年的数据是0.474，为2005年以来最低水平，而自2008年起，基尼系数也在逐年下降。而此前西南财大调查数据显示，中国的2012年的基尼系数为0.61，但无论是民间统计的数据还是官方统计的数据，结果都遭到学术界质疑，仍具有争议性。 本文将根据网络上国家统计局的数据，利用上面给出的公式来计算我国从2002年以来的城镇居民基尼系数，并将计算出的数据与现有数据进行比较。 全球基尼系数 1．考虑两个线性方程组，其系数矩阵如下 1211 11...23211111...1212341,121 1111...3452..............................121111... 12 21n n A A n n n n n ? ???? ?-??????--?? +?? ????==--??????+???? ??-?? ? ?????++-?? 问题的真解均取为[1,1,1,1,...1]T x =，线性方程组的右端项用这个真解计算出来。相应的问题分别称为问题I 和问题II 。请进行如下数值实验： （1） 对问题I 分别用Gauss 消元法，Cholesky 方法，修改的LDLT 算法，追赶法四 种方法求解，其中n=100； （2） 对问题II 分别用Gauss 消去法，列主元Gauss 消去法，不做行交换的列主元 Gauss 消去法求解，其中n=6； （3） 不断增加问题II 的矩阵阶数n=6,8,10,…,20，重复（2）的工作，看看会有什么 问题发生？解释其原因。 （1） Gauss ： 计算程序： n=100; A=2*eye(n); for i=1:n-1 A(i+1,i)=-1; A(i,i+1)=-1; end b=0; b(1)=1; b(100)=1; [x,XA]=GaussJordanXQ(A,b); Gauss 消元法源程序： %用Gauss 消元法解线性方程组 function [x,XA]=GaussJordanXQ(A,b) N = size(A); n = N(1); for i=1:(n-1) for j=(i+1):n if(A(i,i)==0) disp('对角元素为0！'); %防止对角元素为0 return; end l = A(j,i); m = A(i,i); A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; %消元方程 b(j)=b(j)-l*b(i)/m; end end x=SolveUpTriangle(A,b); %通用的求上三角系数矩阵线性方程组的函数XA = A; %消元后的系数矩阵 （SolveUpTriangle.m）解上三角方程组源程序：%解上三角方程组 function x=SolveUpTriangle(A,b) N = size(A); n = N(1); x(n)=b(n)/A(n,n); for k=n-1:1 s=0; for i=k+1:n s=s+A(k,i)*x(i); end x(k)=(b(k)-s)/A(k,k); end 结果： x1=[0,0,0,…..0,0,1]T x2=[0,0,0,…..0,0,0.3820]T x3=[0,0,0,…..0,0,0.9900]T x4=[1,1,1,…..1,1,1]T Cholesky： 《计算方法》课程平时作业 （2010-2011学年第一学期） 作业(考试前交, 给出证明或计算过程、计算程序及计算结果) 1. 对向量()1 2T n x x x x = 定义 1211 , max , n k k k n k x x x x x ∞ ≤≤==== ∑ 设A 是n n ?矩阵，规定 1111 max x A Ax ==，1max x A Ax ∞ ∞∞ ==，2221 max x A Ax == 证明 11111 2max (),max (), . n n kj jk j n j n k k T A a A a A A A λ∞ ≤≤≤≤=====∑∑列范数行范数是最大特征值 2. 用简单迭代法（即不动点迭代法）求方程 32210200x x x ++-= 在1x =附近的根. 要求给出不动点方程、程序和运行结果. 3. 用Newton 迭代法求方程 32210200x x x ++-= 的一个正根，计算结果精确到7位有效数字. 要求给出程序和运行结果. 4. 用牛顿迭代法求方程310x x --=在01x =附近的根. 要求给出程序和运行结果. 5. 证明迭代格式 ()212 33k k k k x x a x x a ++= +， 00,0a x >> . 6. 编写用全主元Gauss 消去法解线性方程组的程序，并求解 12345123451234512345 123450.024******* 42334332416 34418 x x x x x x x x x x x x x x x x x x x x x x x x x -+-+=??-++++=?? +++-=??-++++=??+-++=? 7. 用追赶法解线性方程组 12345210001121000012100001210000120x x x x x -????????????--????????????=--??????--??????????? ?-?????? 要求给出程序和运行结果. 8．给定线性方程组 1212 21 32x x x x +=-?? +=? 问用雅可比迭代法和Gauss-Seidel 迭代法求解是否收敛？ 9. 设有线性方程组 123521121422023103x x x -?????? ??????-=??????????? ?-?????? （1）考察用Jacobi 迭代法和Gauss-Seidel 迭代法解此方程组的收敛性； （2）分别用Jacobi 迭代法和Gauss-Seidel 迭代法解此方程组，要求当 (1)()410k k x x +--<时迭代终止. 给出求解程序和迭代次数及结果. 10．编写幂法程序求矩阵 422251216A ?? ?= ? ??? 按模最大的特征值1λ和对应的特征向量1x 。10(10)ε-= 11. 编写用原点位移加速反幂法程序求矩阵 ?? ?? ? ?????= 3 1- 2-1- 4 3 2- 3 7 A第六章计算机的运算方法(含答案)
数值计算方法大作业
计算方法上机作业
西工大计算方法作业答案
2020年奥鹏吉大网络教育《计算方法》大作业解答
第六章 计算方法简介
西安交通大学计算方法B大作业资料
计算方法上机实习题大作业(实验报告).
计算方法作业2
计算方法大作业非线性方程求根的新方法
工程计算方法及软件应用--本科生考查大作业
第六章数值计算方法举例
北京科技大学计算方法大作业
西交计算方法A上机大作业
计算方法上机作业集合
数值计算方法大作业
计算方法作业第六章
计算方法平时作业
- 第六章货币资金的业务流程与核算方法设计.pptx
- 西北工业大学计算方法第六章作业答案
- 计算方法第六章print
- 6成本会计学-讲义第六章：产品成本计算方法——基本方法
- 数值计算方法-复习-第六章
- 第六章 容器计量的计算方法
- 计算方法第六章迭代法
- 计算方法第六章(迭代法)分析
- 计算方法第六章作业【精选】
- 计算方法第六章
- 计算方法第六章数值积分与数值微分
- 计算方法第六章作业
- 计算方法引论-第六章
- 第六章AHP决策分析的基本原理与计算方法2.pptx
- 计算方法第六章解线性方程组的迭代法
- 计算机组成原理第六章计算机的运算方法
- 第六章计算方法
- 计算方法作业第六章
- 计算方法第六章(迭代法)
- 第六章 计算方法简介