计算方法上机答案
上海电力学院
数值分析上机实验报告
题目：数值分析上机实验报告
学生姓名：11111111111
学号：111111*********
专业：1111
2013年12月30日
数值计算方法上机实习题
1． 设?+=1
05dx x
x I n
n ， （1） 由递推公式n I I n n 1
51+
-=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用n
I I n n 51
511+-=-，计算0I ；
（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
(1) 解答：n=0，
0.1823)05ln()15ln()5(5151510101
0=+-+=++=+=+=???x d x
dx x dx x x I n
n
这里可以用for 循环，while 循环，根据个人喜好与习惯：
for 循环程序： While 循环程序： I=0.1823; I=0.1823; for n=1:20 i=1;
I=(-5)*I+1/n; while i<21 End I=(-5)*I+1/i; I i=i+1; fprintf('I20=%f',I) end I = -2.0558e+009 >> I
I20=-2055816073.851284>> I = -2.0558e+009 (2) 粗略估计I 20： Mathcad 计算结果： for 循环程序： While 循环程序： >> I=0.007998; I=0.007998; >> for n=1:20 n=1;
I=(-0.2)*I+1/(5*n); while n<21
End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> I
I =0.0083
(3) 算法误差分析：
计算在递推过程中传递截断误差和舍入误差 第一种算法：（从1——>20）
*
000
e I I
=-
*
**21111120
11
5(5)5()555n n n n n n n n n n e I I I I I I e e e n n
------=-=-+
--+=-===
1
x x 205x +?????
d 7.99810
3
-?=
误差放大了5n 倍，算法稳定性很不好； 第二种算法：（从20——>1）
*
n n n
e I I =-
*
**111111111()()555555n n n n n n n
n e I I I I I I e n n ---=-=-+--+=-=
0111...()55
n
n
e e e ===
误差在逐步缩小，算法趋近稳定，收敛。
2． 求方程0210=-+x e x 的近似根，要求4
1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法；
function [t i]=erfenfa(a,b) f=@(x)( exp(x)+10*x-2) t=(a+b)./2; i=0;
while abs(f(t))>0.001 if f(a)*f(t)<0 b=t;t=(a+b)/2; elseif f(b)*f(t)<0 a=t;t=(b+a)/2; end i=i+1; end 结果： t = 0.0906 i = 11
（2） 取初值00=x ，并用迭代10
21
x k e x -=+；
function x=diedai(x0) %x0初值 x=x0;
for i=1:10000
y=(2-exp(x))./10;x=y;y=(2-exp(x))./10; if abs(x-y)<5*10^(-4) disp('迭代次数'); 2*i break; end end
结果：
ans =
6 x =
0.090639135859584
（3） 加速迭代的结果(艾特肯Aitken 加速方法)；
function [y m]=aitken(a) func=@(x)(( 2-exp(x))./10) x(1)=a; wucha=1;m=1;
while wucha> 5*10^(-4) p(m+1)=func(x(m));
q(m+1)=func(p(m+1));
x(m+1)=q(m+1)-((q(m+1)-p(m+1))^2)./(q(m+1)-2*p(m+1)+x(m));
wucha=abs(x(m+1)-x(m)); m=m+1; if m>1000 break; end end
format long y=x(m-1); m=m-1;
运行结果y =0.090483741803596
m =2
（4） 取初值00 x ，并用牛顿迭代法；
function x=newtondiedai(x0) x=x0; for i=1:100
y=x-(exp(x)+10*x-2)./(exp(x)+10);x=y; y=x-(exp(x)+10*x-2)./(exp(x)+10); if abs(x-y)<0.0001 disp('迭代次数'); i break; end end 运行结果 迭代次数
i = 2 x =
0.0905
（5） 分析绝对误差。
通过指令求得方程精确解的近似解：
>> solve('exp(x)+10*x-2=0') ans =
-lambertw(1/10*exp(1/5))+1/5 >> format long
>> -lambertw(1/10*exp(1/5))+1/5 ans =
0.09052510130725 小结：
所以我们可以看到，在要求同一运算精度的情况下，采用二分法运算了11次，采用题设的迭代方法运算了6次，采用加速迭代法只运算了2次，采用牛顿迭代法需运算2次。也就是说加速迭代速度都是超线性收敛的，而简单迭代法是线性收敛的。而二分法收敛速度较慢,所以在工程上不经常使用。
3．钢水包使用次数多以后，钢包的容积增大，数据如下： x 2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
y
6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76
试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型） 解答：因为不知道其函数形式，可先plot 而后确定使用哪种模型比较合适。 函数图形程序：
x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76]; plot(x,y ，‘b*’)
与指数函数趋势类似（但是趋势不同，一个趋于递减，一个趋于递增，使用倒数），故拟合模型为： b
x
y ae
=
两边同时取对数：
ln ln b
y a
x
=+
用此模型进行数据拟合：
1,ln n n y x ?? ???
x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76]; t=1./x;
s=polyfit(t,log(y),1) s =
-1.1107 2.4578
即是： 2.457811.679a e == 1.1107b =- 最终拟合曲线为：
1.1107
11.679x
y e
-=
将此拟合曲线和源数据进行对比： 主程序：
x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76]; plot(x,y,'ro') hold on
lims1=[2,16];
fplot('11.679*exp(-1.1107/x)',lims1) hold off
均方差的求解
x=[2:16];
y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76]; f(x)=11.679*exp(-1.1107./x); c=0;
for i=1:15 a=y(i); b=x(i);
c=c+(a-f(b))^2;
end
averge=c/15
averge = sqrt(averge) 0.2438
4．设??
?
?
??
?
???
?
?----------------=41010014101001410110141001014100
1014A ，?????????? ??--=625250b ，b x =A 分析下列迭代法的收敛性，并求42
110-+≤-k
k x x 的近似解及相应的迭代次数。
1） JACOBI 迭代；
2） GAUSS-SEIDEL 迭代； 3） SOR 迭代（95.0,
95.1,334.1=ω）
。 解答：
（1） JACOBI 迭代； 定义函数。
function tx=jacobi(A,b,imax,x0,tol) %jacobi 迭代 %初始值x0,次数imax,精度tol del=10^-10;
tx=[x0];n=length(x0); for i=1:n
dg=A(i,i); if abs(dg) disp('对角元太小'); return end end for k=1:imax %jacobi 循环 for i=1:n sm=b(i); for j=1:n if j~=i sm=sm-A(i,j)*x0(j); end end x(i)=sm/A(i,i); end tx=[tx;x]; if norm(x-x0) x0=x; end end 保存m文件。 主窗口输入以下程序段： >>A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]; x0=[0 0 0 0 0 0]; imax=100; tol=10^-4; tx=jacobi(A,b,imax,x0,tol); for j=1:size(tx,1) fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)); end 运行结果： 1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2 0.000000 1.250000 -0.500000 1.250000 -0.500000 1.500000 3 0.625000 1.000000 0.500000 1.000000 0.500000 1.250000 4 0.500000 1.656250 0.312500 1.656250 0.312500 1.750000 5 0.828125 1.531250 0.765625 1.531250 0.765625 1.656250 6 0.765625 1.839844 0.679688 1.839844 0.679688 1.882813 7 0.919922 1.781250 0.890625 1.781250 0.890625 1.839844 8 0.890625 1.925293 0.850586 1.925293 0.850586 1.945313 9 0.962646 1.897949 0.948975 1.897949 0.948975 1.925293 10 0.948975 1.965149 0.930298 1.965149 0.930298 1.974487 11 0.982574 1.952393 0.976196 1.952393 0.976196 1.965149 12 0.976196 1.983742 0.967484 1.983742 0.967484 1.988098 13 0.991871 1.977791 0.988895 1.977791 0.988895 1.983742 14 0.988895 1.992415 0.984831 1.992415 0.984831 1.994448 15 0.996208 1.989639 0.994820 1.989639 0.994820 1.992415 16 0.994820 1.996462 0.992923 1.996462 0.992923 1.997410 17 0.998231 1.995167 0.997583 1.995167 0.997583 1.996462 18 0.997583 1.998349 0.996699 1.998349 0.996699 1.998792 19 0.999175 1.997745 0.998873 1.997745 0.998873 1.998349 20 0.998873 1.999230 0.998460 1.999230 0.998460 1.999436 21 0.999615 1.998948 0.999474 1.998948 0.999474 1.999230 22 0.999474 1.999641 0.999282 1.999641 0.999282 1.999737 23 0.999820 1.999509 0.999755 1.999509 0.999755 1.999641 24 0.999755 1.999832 0.999665 1.999832 0.999665 1.999877 25 0.999916 1.999771 0.999886 1.999771 0.999886 1.999832 26 0.999886 1.999922 0.999844 1.999922 0.999844 1.999943 27 0.999961 1.999893 0.999947 1.999893 0.999947 1.999922 28 0.999947 1.999964 0.999927 1.999964 0.999927 1.999973 29 0.999982 1.999950 0.999975 1.999950 0.999975 1.999964 （2）GAUSS-SEIDEL迭代； function tx=gseidel(A,b,imax,x0,tol) %gseidel迭代 %初始值x0,次数imax,精度tol del=10^-10; tx=[x0];n=length(x0); for i=1:n dg=A(i,i); if abs(dg) disp('对角元太小'); return end end for k=1:imax %G-S循环 x=x0; for i=1:n sm=b(i); for j=1:n if j~=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); end tx=[tx;x]; if norm(x-x0) return else x0=x; end end 主程序： A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]; x0=[0 0 0 0 0 0]; imax=100; tol=10^-4; tx=gseidel(A,b,imax,x0,tol); for j=1:size(tx,1) fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5 ),tx(j,6)) end ω）。 （3）SOR迭代（95 .1 334 = 95 .0 , .1 , function tx=sor(A,b,imax,x0,tol,w) %sor迭代 %初始值x0,次数imax,精度tol del=10^-10; tx=[x0];n=length(x0); for i=1:n dg=A(i,i); if abs(dg) disp('对角元太小'); return end end for k=1:imax %G-S循环 x=x0; for i=1:n sm=b(i); for j=1:n if j~=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); x(i)=w*x(i)+(1-w)*x0(i); end tx=[tx;x]; if norm(x-x0) return else x0=x; end end 主程序： A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]; x0=[0 0 0 0 0 0]; imax=100; tol=10^-4; w=1.334; tx=sor(A,b,imax,x0,tol,w); for j=1:size(tx,1) fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end 分析： 从题中可以看出，系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛；又因系数矩阵对称正定，且0＜w＜2，则SOR迭代法也收敛。 三种迭代法的结果分析及比较：（1）Jacobi法，其设计思想是将线性方程组的求解归结 为对角方程组求解过程的重复，计算规则简单而易于编写汁算程序。但，我们通常要设置最大迭代次数N ，以防止迭代过程不收敛或者收敛速度过于缓慢。（2）Gauss-Seidel 充分利用新信息进行计算，其迭代的效果比Jacobi 迭代好。（3）SOR 迭代法，①计算公式简单、编制程序容易，是求解大型稀疏方程组的一种有效方法，且超松弛法收敛速度比Gauss-Seidel 迭代和Jacobi 迭代都要快。② 从对松弛因子三种不同选择的分析发现，若松弛因子选择得当，能显著地提高收敛速度。 5．用逆幂迭代法求??? ? ? ??=111123136A 最接近于11的特征值和特征向量，准确到310-。 解答： 调用函数： function [mt,my]=maxtr(A,p,ep) n=length(A); B=A-p*eye(n); v0=ones(n,1); k=1; v=B*v0; while abs(norm(v,inf)-norm(v0,inf))>ep %norm(v-v0)>ep k=k+1; q=v; u=v/norm(v,inf) v=B*u; v0=q; end mt=1/norm(v,inf)+p my=u 主程序： A=[6 3 1;3 2 1;1 1 1]; maxtr(A,11,0.001) 6．用经典R-K 方法求解初值问题 （1）???-+-='++-='x x y y y x y y y sin 2cos 22sin 22212 211 ，]10,0[∈x ， ???==3)0(2)0(21y y ； （2）???-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212211，]10,0[∈x ， ???==3 )0(2 )0(21y y 。 和精确解???+=+=--x e x y x e x y x x cos 2)(sin 2)(21比较，分析结论。 解答： （1）主程序： function ydot=lorenzeq(x,y) ydot=[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)]; 主窗口输入： x=0:10; y1=2*exp(-x)+sin(x); y2=2*exp(-x)+cos(x); [x,y]=ode45('lorenzeq',[0:10],[2;3]); fprintf(' x y(1) y1 y(2) y2\n') for j=1:length(y) fprintf('%4d %.4f %.4f %.4f %.4f\n',x(j),y(j,1),y1(j),y(j,2),y2(j)) end 结果： x y(1) y1 y(2) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5775 1.577 2 1.2758 1.2761 2 1.1802 1.1800 -0.1457 -0.1455 3 0.2406 0.2407 -0.8903 -0.8904 4 -0.7202 -0.7202 -0.6170 -0.6170 5 -0.9454 -0.9454 0.2971 0.2971 6 -0.2745 -0.2745 0.9652 0.9651 7 0.6589 0.6588 0.7557 0.7557 8 0.9901 0.9900 -0.1449 -0.1448 9 0.4124 0.4124 -0.9109 -0.9109 10 -0.5440 -0.5439 -0.8389 -0.8390 （2）主程序： function ydot=lorenzeq1(x,y) ydot=[-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)]; 主窗口输入： x=0:10; y1=2*exp(-x)+sin(x); y2=2*exp(-x)+cos(x); [x,y]=ode45('lorenzeq1',[0:10],[2;3]); fprintf(' x y(1) y1 y(2) y2\n') for j=1:length(y) fprintf('%4d %.4f %.4f %.4f %.4f\n',x(j),y(j,1),y1(j),y(j,2),y2(j)) end 运行结果： x y(1) y1 y(2) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.577 2 1.5772 1.2759 1.2761 2 1.1800 1.1800 -0.1455 -0.1455 3 0.2407 0.2407 -0.890 4 -0.8904 4 -0.7202 -0.7202 -0.6169 -0.6170 5 -0.9454 -0.9454 0.2972 0.2971 6 -0.2745 -0.2745 0.9648 0.9651 7 0.6588 0.6588 0.7554 0.7557 8 0.9900 0.9900 -0.1448 -0.1448 9 0.4124 0.4124 -0.9106 -0.9109 10 -0.5439 -0.5439 -0.8389 -0.8390 小结： 四阶RungeKutta 方法得到的结果已很接近精确解，证明这种迭代方法精确度很好，是一种有效的算法。 7．用有限差分法求解边值问题（h=0.1）： ?? ?==-=+-''1 )1()1(0 )1(2y y y x y . 解答： 主程序： h=0.1; n=(1-(-1))/h+1; x(1)=-1;x(n-1)=1; y(1)=1;y(n-1)=1; for i=1:n-1 x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)^2); B(i)=2/(h^2)+q(i); end for i=1:n-2 C(i)=-1/(h^2); end H=diag(B)+diag(C,1)+diag(C,-1); g(1)=0+1/(h^2); g(n-1)=0+1/(h^2); for i=2:n-2 g(i)=0; end y=H\g' 运行结果： y = 0.9027 0.8235 0.7592 0.7074 0.6661 0.6338 0.6095 0.5922 0.5814 0.5767 0.5778 0.5846 0.5974 0.6163 0.6420 0.6752 0.7167 0.7680 0.8308 0.9072 小结： 有限差分法是将微分方程离散化为差分方程进行求解，而求解差分方程就是解一个三对角矩阵线性方程组，从以上程序可以看出编程很容易实现用追赶法求解三对角矩阵线性方程组，且具有较高的精度。 8．用函数bx a y sin =拟合数据 x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 y 0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.5 解：此为一个非线性最小二乘问题，按照最小二乘原理，应选取参数b a ,使得表达式 ∑==8 1i 2 i i )y -(asinbx I 达到极小值。用matlab 仿真需要借助fminsearch 函数。 完整的matlab 程序如下： 解答：非线性最小二乘问题： x=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8]; y=[0.6 1.1 1.6 1.8 2 1.9 1.7 1.4]; 第一步：创建 fitfun 函数 fitfun.m function err=fitfun(c,x,y) a=c(1); b=c(2); err=y-a*sin(b*x); err=err'*err; 第二步：建立nlfit.m,利用fminsearch 寻找fitfun 极小值： function [err,a,b]=nlfit(x,y) if nargin<2 %判断输入变量个数 x=[1:8]'/10; %输入数据 y=[0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3]'; end c=fminsearch(@fitfun,[0;0],optimset,x,y); %a,b 都从0开始搜索 fprintf('The nonlinear least square fitting y=a*sin(b*x) for data\n\n'); fprintf('%6.1f',x); fprintf('\n'); %换行 fprintf('%6.1f',y); %6位浮点数表示保留1位小数 fprintf('\n\n is\n\t y=%7.4f*sin(%7.4f*x)\n\n',c(1),c(2)); z=linspace(x(1),x(end),100); %x(1)至x(end)分隔100份,为了画出曲线 plot(x,y,'r+',z,c(1)*sin(c(2)*z),'b-.') 结果： >> nlfit The nonlinear least square fitting y=a*sin(b*x) for data 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3 is y= 1.9751*sin( 3.0250*x) 最终输出图如图所示： 0.1 0.20.30.40.50.60.70.8 0.51 1.5 2 9.拟合形如 ()1a bx f x cx +≈ +的函数的一种快速方法是将最小二乘法用于下列问题： ()(1)f x cx a bx +≈+，试用这一方法拟合下面的表格给出的中国人口数据。 中国人口数据 次序 年份 人口（亿） 第一次 1953 5.82 第二次 1964 6.95 第三次 1982 10.08 第四次 1990 11.34 第五次 2000 12.66 x=[1953 1964 1982 1990 2000]'; >> y=[5.82 6.95 10.08 11.34 12.66]'; >> x1=ones(5,1); >> A=[x1 x x.*y]; >> z=A\y z = 2.9456 -0.0014 0.0005 因此，方程为：y=( 2.9456-0.0014*x)/(1+0.0005*x) 10、已知20世纪美国人口的统计数字如表所示（单位：百万） 表美国人口统计数据 年份1900 1910 1920 1930 1940 1950 1960 1970 1980 1990人口76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 解:由logisitc模型， p M p r t P) 1( )(- = ' 。其中r为固有增长率，M为最大人口容 量。 解得 rt e p M M t P - - + = )1 ( 1 )( 0。 由程序： t=1900:10:1990; p=[76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4]; n=10; dp=p(2:n)-p(1:n-1); dt=t(2:n)-t(1:n-1); dpdt=(dp./dt)./p(2:n); f=polyfit(p(2:n),dpdt,1); disp('求人口增长率r：'),r=f(2) disp('求最大人口容量M：'),s=-f(1);M=r/s p1=M./(1+(M/p(1)-1).*exp(-r.*(t-1900))); plot(t,p,'*',t,p1) 运行结果为： 求人口增长率r： r = 0.0164 求最大人口容量M： M =659.8734 美国人口增长模型 综上所诉：本题利用logisitc模型拟合，先拟合 ) 1( M p r 。 人口的长期预测是一件复杂而又困难的事情,但中短期的预报却是现实可能的,。人口增长都不是线性关系,美国人口在局部时段是指数增长的, 但长期趋势却是Logistic 过程。 . / 《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：XXX 二〇一八年五月二十五日 前言 通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。 以下为本次上机实验报告，按照实验内容共分为六部分。 实验一： 一、实验名称及题目： Newton 迭代法 例2.7(P38):应用Newton 迭代法求在附近的数 值解，并使其满足. 二、解题思路： 设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交 点的横坐标) (') (0001x f x f x x - =,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标) (') (1112x f x f x x - =称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把) (') (1n n n n x f x f x x - =+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。 三、Matlab 程序代码： function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0); y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1; while abs(x1-x0)>=tol x0=x1; y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; end x=double(x1) K 四、运行结果： 数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名： 课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m ) 则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。 2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量 （1）二分法 （局部，大图不太看得清，故后面两小题都用局部截图） （2）迭代法 （3）牛顿法 顺序消元法 #include for(k=p+1;k 计算方法上机报告 姓名： 学号： 班级：能动上课班级： 题目及求解： 一、对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=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 算法结构 ;0=s ?? ? ??+-+-+-+= 681581482184161n n n n t n ； for 0,1,2,,n i =??? if 10m t -≤ end; 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=7，此时求解得： s =3.1415926536； 若保留30位有效数字时，则n=22, 此时求解得： s =3.8。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 二、某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示： 第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1，其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________') 三、结果 从结果可以看出有效位数是6位。 感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。 数值计算方法上机实习题 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+ -=-，从I 0=0.1824, 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n I I n n 51 5111+- =--，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 答：第一个算法可得出 e 0=|I 0?I 0 ?| e n =|I n ?I n ?|=5n |e 0| 易知第一个算法每一步计算都把误差放大了5倍，n 次计算后更是放大了5n 倍，可靠性低。 第二个算法可得出 e n =|I n ?I n ?| e 0=(15 )n |e n | 可以看出第二个算法每一步计算就把误差缩小5倍，n 次后缩小了5n 倍，可靠性高。 2． 求方程0210=-+x e x 的近似根，要求41105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 计算根与步数程序： fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1); fprintf('root=%6.8f ,n=%d \n',root,n); 计算结果显示： root=0.09057617 ,n=11 （2） 取初值00=x ，并用迭代10 21 x k e x -=+； （3） 加速迭代的结果; （4） 取初值00 x ，并用牛顿迭代法； 《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：XXX 二〇一八年五月二十五日 前言 通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。 以下为本次上机实验报告，按照实验内容共分为六部分。 实验一： 一、实验名称及题目： Newton 迭代法 例2.7(P38):应用Newton 迭代法求 在 附近的数值解 ，并使其满足 . 二、解题思路： 设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标) (') (0001x f x f x x - =,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标) (') (1112x f x f x x - =称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把 ) (') (1n n n n x f x f x x - =+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。 三、Matlab 程序代码： function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0); y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1; while abs(x1-x0)>=tol x0=x1; y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; end x=double(x1) K 四、运行结果： 实验二： 数值计算方法I 上机实验考试题（两题任选一题） １．小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. A. 建立火箭升空过程的数学模型（微分方程）； B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． ２．小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为 0)0(,0,01222==≤≤-+?? ? ??-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻． 今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））： 现有两种估计比例系数k 的方法： 1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。 2．用这组数据拟合一个k ． 请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）． 《数值计算方法》上机实验报告华北电力大学 实验名称数值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 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; 计算方法 一、填空题 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 算法语言的基本符号分为 字母 、 数字 、 逻辑值、 定义符四大 实习题1 1用两种不容的顺序计算644834.11000 12≈∑=-n n ，分析误差的变化 （1）顺序计算 源代码： #include if(n%1000==0) printf("sum[%d]=%-30f",n,sum); if(n<=1)break; n--; } printf("sum[%d]=%f\n",n,sum); } 结果： 2已知连分数 )) / /(... /( 3 2 2 1 1 n n b a a b a b a b f + + + + = 利用下面的方法计算f: 1 1)0 ,..., 2 ,1 ( , d f n n i d a b d b d i i i i n n = - - = + = = + + 写一个程序，读入n, n n b a,,计算并打印f 源代码： #include 本科实验报告 课程名称：计算机数值方法 实验项目：方程求根、线性方程组的直接解 法、线性方程组的迭代解法、代数插值和最 小二乘拟合多项式 实验地点：行勉楼 专业班级： ******** 学号： ********* 学生姓名： ******** 指导教师：李誌，崔冬华 2016年 4 月 8 日 y = x*x*x + 4 * x*x - 10; return y; } float Calculate(float a,float b) { c = (a + b) / 2; n++; if (GetY(c) == 0 || ((b - a) / 2) < 0.000005) { cout << c <<"为方程的解"<< endl; return 0; } if (GetY(a)*GetY(c) < 0) { return Calculate(a,c); } if (GetY(c)*GetY(b)< 0) { return Calculate(c,b); } } }; int main() { cout << "方程组为：f(x)=x^3+4x^2-10=0" << endl; float a, b; Text text; text.Getab(); a = text.a; b = text.b; text.Calculate(a, b); return 0; } 2.割线法： // 方程求根（割线法）.cpp : 定义控制台应用程序的入口点。// #include "stdafx.h" #include"iostream" 心得体会 使用不同的方法，可以不同程度的求得方程的解，通过二分法计算的程序实现更加了解二分法的特点，二分法过程简单，程序容易实现，但该方法收敛比较慢一般用于求根的初始近似值，不同的方法速度不同。面对一个复杂的问题，要学会简化处理步骤，分步骤一点一点的循序处理，只有这样，才能高效的解决一个复杂问题。 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。 Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear %%%给定求解区间 b=1.5; a=0; %%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f); 计算方法实验报告 班级： 学号： 姓名： 成绩： 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．理解数值计算稳定性的概念． 2．了解数值计算方法的必要性． 3．体会数值计算的收敛性与收敛速度． 三、实验内容 计算dx x x I n n ? += 1 10 ，1,2,,10n = ． 四、算法描述 由 dx x x I n n ? += 1 10 ，知 dx x x I n n ?+=--101110，则 上海电力学院 数值分析上机实验报告 题目：数值分析上机实验报告 学生姓名：11111111111 学号：111111********* 专业：1111 2013年12月30日 数值计算方法上机实习题 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+ -=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用n I I n n 51 511+-=-，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 (1) 解答：n=0， 0.1823)05ln()15ln()5(5151510101 0=+-+=++=+=+=???x d x dx x dx x x I n n 这里可以用for 循环，while 循环，根据个人喜好与习惯： for 循环程序： While 循环程序： I=0.1823; I=0.1823; for n=1:20 i=1; I=(-5)*I+1/n; while i<21 End I=(-5)*I+1/i; I i=i+1; fprintf('I20=%f',I) end I = -2.0558e+009 >> I I20=-2055816073.851284>> I = -2.0558e+009 (2) 粗略估计I 20： Mathcad 计算结果： for 循环程序： While 循环程序： >> I=0.007998; I=0.007998; >> for n=1:20 n=1; I=(-0.2)*I+1/(5*n); while n<21 End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> I I =0.0083 (3) 算法误差分析： 计算在递推过程中传递截断误差和舍入误差 第一种算法：（从1——>20） * 000 e I I =- * **21111120 11 5(5)5()555n n n n n n n n n n e I I I I I I e e e n n ------=-=-+ --+=-=== 1 x x 205x +????? d 7.99810 3 -?= 计算方法B 上机报告 第1题 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示： (1)请用合适的曲线拟合所测数据点； (2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图； 问题分析和算法思想： 本题的主要目的是对21个测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值，可以用的插值方法很多：多项式插值、Lagrange 插值、Newton 插值、三次样条插值等。由于数值点较多时，采用高次多项式插值将产生很大的误差，用拉格朗日插值多项式会出现龙格现象。故为了将所有的数据点都用上，且题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值的方法较为合适。 计算光缆长度近似值，只需将每两点之间的距离算出，然后依次相加，所得的折线长度，即为光缆长度的近似值。 光缆长度计算公式： 19 1 k k k l +===∑? ? ? 算法结构： 三次样条算法结构见《计算方法教程》P110。 源程序： clear;clc; x=0: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.90 7.95 8.86 9.81 10.80 10.93]; d=y; plot(x,y,'k.','markersize',15) hold on %%%计算二阶差商 for k=1:2 for i=21:-1:(k+1) d(i)=(d(i)-d(i-1))/(x(i)-x(i-k)); end end %%%假定d的边界条件，采用自然三次样条 for i=2:20 d(i)=6*d(i+1); end d(1)=0; d(21)=0; %%%追赶法求解带状矩阵的m值 a=0.5*ones(1,21); b=2*ones(1,21); c=0.5*ones(1,21); a(1)=0;c(21)=0; u=ones(1,21); u(1)=b(1); r=c; yy(1)=d(1); %%%追的过程 for k=2:21 l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); yy(k)=d(k)-l(k)*yy(k-1); end %%%赶的过程 m(21)=yy(21)/u(21); for k=20:-1:1 m(k)=(yy(k)-r(k)*m(k+1))/u(k); end %%%利用插值点画出拟合曲线 k=1; nn=100; xx=linspace(0,20,nn); l=0; for j=1:nn for i=2:20 if xx(j)<=x(i) k=i; 第一次&第二次上机作业 上机作业： 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); 实验一：误差传播与算法稳定性 实验目的：体会稳定性在选择算法中的地位。 实验内容：考虑一个简单的由积分定义的序列 1 0I ,0,1,10n n x dx n a x ==+? 其中a 为参数，分别对0.05a =及15a =按下列两种方法计算。 方案1：用递推公式 11,1,2,,10n n I aI n n -=-+= 递推初值可由积分直接得 01 ln a I a += 方案2：用递推公式 111 (), ,1,,1n n I I n N N a n -=-+=- 根据估计式 当1n a n ≥ +时，11 (1)(1)(1) n I a n a n <<+++或 当01n a n ≤< +时， 11 (1)(1)n I a n n <≤++ 取递推初值 当1n a n ≥ +时， 11121()2(1)(1)(1)2(1)(1)N N a I I a N a N a a N +≈+=+++++ 当01n a n ≤< +时，111()2(1)(1)N N I I a N N ≈+++ 实验要求：列出结果，并对其稳定性进行分析比较，说明原因。 实验二：非线性方程数值解法 实验目的：探讨不同方法的计算效果和各自特点 实验内容：应用算法（1）牛顿法；（2）割线法 实验要求： （1）用上述各种方法，分别计算下面的两个例子。在达到精度相同的前提下，比较其迭代次数。 （I ）3 1080x x +-=，取00x =； (II) 2 2 81(0.1)sin 1.060x x x -+++=，取00x =； （2） 取其它的初值0x ，结果如何？反复选取不同的初值，比较其结果； （3） 总结归纳你的实验结果，试说明各种方法的特点。 实验三：选主元高斯消去法----主元的选取与算法的稳定性计算方法上机实验报告
数值分析上机作业
计算方法上机题答案
西安交通大学计算方法B上机报告
数值分析上机题目详解
数值计算方法上机实习题
计算方法上机实验报告
数值计算方法I上机实验考试题
《数值计算方法》上机实验报告
(完整版)数值计算方法上机实习题答案
计算方法试题库讲解
计算方法与实习上机题答案
太原理工大学数值计算方法实验报告
(完整版)哈工大-数值分析上机实验报告
计算方法上机实习题大作业(实验报告).
计算方法实验报告格式
计算方法上机答案
计算方法B上机报告
计算方法上机作业集合
计算方法上机实验题1~6