计算方法B上机实习报告
计算方法B 上机实习报告
1.对以下和式计算：01421181
84858616n n S n n n n ∞
=??
=--- ?++++??∑
，要求： (1)若只需保留11个有效数字，该如何进行计算；
(2)若要保留30个有效数字，则又将如何进行计算；
问题分析：
在该题中S 的每一项14211()1681848586
n n
s n n n n =
---++++存在两个相近的数相减的问题，因此为了避免有效数字损失，最好是改变运算顺序，分别将正数和负数相加，然后再将其和相加。另外，sn 中有多个负数相加，可以按照绝对值递增的顺序求和，以减少舍入误差的影响。同时，为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需 要的迭代次数，然后采用倒序相加的方法实现。
程序实现：
clear;clc;
m=input('请输入要保留的有效数字位数：'); s1=0; s2=0; k=0; s=1;
%%%%判断多需要的迭代次数 while s>=0.5*10^-(m-1)
s=4/(16^k*(8*k+1))-(2/(16^k*(8*k+4))+1/(16^k*(8*k+5))+1/(16^k*(8*k+6))); k=k+1; end
%%%%正负数分别按照绝对值递增的顺序倒序相加 for n=(k-1):-1:0 a1=4/(16^n*(8*n+1)); a2=2/(16^n*(8*n+4)); a3=1/(16^n*(8*n+5)); a4=1/(16^n*(8*n+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=s1-s2; S=vpa(S,m)
运算结果：
总结心得：
在计算求和问题中，应特别注意相近数相减的问题，这样会造成有效数字灾难性的损失。另外在两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，因此要按照从小到大的顺序相加。在上题计算中分别对正负相采用倒序相加，这样就有效的避免了“大数吃小数”的问题。
2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)
(1)请用合适的曲线拟合所测数据点；
(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；
问题分析：
本题的主要目的是对测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值。由于数值点较多时，使用拉格朗日差值多项式会出现龙格现象。为了将所有的数据点都用上，采用分段差插法，本题使用三次样条插值。
算法思想：
M记样条函数在每个子区间上是三次多项式，它的二阶导数必是一次多项式。若用
i
在i x 处的二阶导数''(x )i s 。则在区间1[x ,x ]i i - 上 ''1
1(x)i i i i i
x x x x s M M
h h ----=+ 式中 1i i i h x x -=- （1） 对上式进行两次积分得
332211
111()()()()()
6666i i i i i i i i i i i i i i i i
x x x x h x x h x x s x M M y M y M h h h h ---------=++-+- （2） 它的一阶导数为
22'
111
1()()()226
i i i i i i i i i i i i x x x x y y M M s x M M h h h h --------=-++- （3）
'()i s x 满足连续性条件，即
''
(x )()i i s s x -+= 1,2,n i =-…
,1 上式和（3）式得
1111
111636i i i i i i i
i i i i i i
h h h h y y y y M M M h h +++--+++--++=- （4） 用差商记号，并记
11i i i i h h h λ++=
+ 1
i
i i i h h h μ+=+
（4）式可以写成
111126[,,]i i i i i i i i M M M y x x x μλ-+-+++= 1,2,1i n =-…
, 方程组可以写成如下形式
0011111
111222
2n n n n n
n n M d M d M d M d λμλμλμ----??????
??? ? ??? ? ??? ?= ??? ? ??? ?
??? ??
?????
116[,,]i i i i d y x x x -+= i 1,2
,n =…, 自然样条插值条件为
000,0,0,0n n d d λμ====
在估计河底光缆长度时使用第一类线积分
20
19
k k
k l ds +====∑???
长度
程序实现：
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;
break;
else
k=i+1;
end
end
h=1;
xbar=x(k)-xx(j);
xmao=xx(j)-x(k-1);
s(j)=(m(k-1)*xbar^3/6+m(k)*xmao^3/6+(y(k-1)-m(k-1)*h^2/6)*xbar+(y(k)-m(k)*h^2/6)*xmao)/h; sp(j)=-m(k-1)*(x(k)-xx(j))^2/(2*h)+m(k)*(xx(j)-x(k-1))^2/(2*h)+(y(k)-y(k-1))/h-(m(k)-m(k-1))*h/6; l(j+1)=(1+sp(j)^2)^0.5*(20/nn)+l(j);%利用第一类线积分求河底光缆的长度
end
%%%绘图
plot(xx,s,'r-','linewidth',1.5)
grid
disp(['所需光缆长度为',num2str(l(nn+1)),'米'])
运行结果：
总结心得：
采用三次样条插值对数据进行拟合时，可以有效避免龙格现象。在本题的计算中采用自然三次样条函数的边界条件。在解线性方程组时使用了追赶法求解带状矩阵，在求解三对角矩阵时追赶法计算速度快，是一种求解线性方程组的有效手段。在估计河底光缆长度时使用第一类线积分。本题计算中间变量非常多，在调试的过程中遇到了一些麻烦，这更加使我认识到在编程的过程中由不得一点儿马虎，在下面问题的处理中我变得更加小心。
3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温
问题分析：
对一天的气温进行数据拟合，可以考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数。问题的难度是求解各种拟合函数的系数。利用法方程求解最小二乘系数时，方程的解不够稳定，本题利用正交化方法。
程序实现：
clear;clc;
x=0:24;
y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16];
m=length(x);
n=input('请输入函数的次数');
plot(x,y,'k.',x,y,'-')
grid;
hold on;
n=n+1;
G=zeros(m,n+1);
G(:,n+1)=y';
c=zeros(1,n);%建立c来存放σ
q=0;
f=0;
b=zeros(1,m);%建立b用来存放β
%%%形成矩阵G
for j=1:n
for i=1:m
G(i,j)=x(1,i)^(j-1);
end
end
%%%建立矩阵Qk
for k=1:n
for i=k:m
c(k)=G(i,k)^2+c(k);
end
c(k)=-sign(G(k,k))*(c(k)^0.5);
w(k)=G(k,k)-c(k);%建立w来存放ω
for j=k+1:m
w(j)=G(j,k);
end
b(k)=c(k)*w(k);
%%%变换矩阵Gk-1到Gk
G(k,k)=c(k);
for j=k+1:n+1
q=0;
for i=k:m
q=w(i)*G(i,j)+q;
end
s=q/b(k);
for i=k:m
G(i,j)=s*w(i)+G(i,j);
end
end
end
%%%求解三角方程Rx=h1
a(n)=G(n,n+1)/G(n,n);
for i=n-1:(-1):1
for j=i+1:n
f=G(i,j)*a(j)+f;
end
a(i)=(G(i,n+1)-f)/G(i,i); %a（i）存放各级系数
f=0;
end
a
%%%拟合后的回代过程
p=zeros(1,m);
for j=1:m
for i=1:n
p(j)=p(j)+a(i)*x(j)^(i-1);
end end
plot(x,p,'r*',x,p,'-'); E2=0;%用E2来存放误差 %%%误差求解 for i=n+1:m
E2=G(i,n+1)^2+E2; end E2=E2^0.5; disp('误差为'); disp(E2); t=0 for i=1:m t=t+p(i); end
t=t/m; %%%平均温度
disp(['平均温度为',num2str(t),'℃'])
运行结果：
二次函数时，系数a1=8.3063,a2=2.6064,a3=-0.0938。误差为16.7433，平均温度为21.2℃。 函数形式为 2
()8.3063 2.60640.0938p x x x =+-
三次函数时，系数a1=13.3880,a2=-0.2273,a3=0.2075,a4=-0.0084,误差为11.4482，平均温度为21.2℃。
函数形式 2
3
()13.38800.22730.20750.0084p x x x x =-+-
四次函数时，系数a1=16.7939,a2=-3.7050,a3=0.8909,a4=-0.0532,a5=0.0009,误差为7.6838，平均温度为21.2℃。
函数形式为 2
3
4
()16.7939
3.7050
0.8909
0.05320.0009
p x x x x x =-+-+
当为指数函数时，假定指数函数的形式为2
123a a x a x y e ++= 。只需将y 值求对数，其主体部分不变，求得的系数为a1=2.3835,a2=0.1251,a3=-0.0044,误差为14.6867，平均温度为20.9399℃。 函数形式为 2
2.38350.12510.0044x x y e +-=
四种函数图像的对比：红色表示指数函数，绿色表示二次函数，青色表示三次函数，紫色表示四次函数，蓝色的折线表示数据点。
总结心得：
同上面的结果可以看出，随着次数的增加误差在逐渐减小，拟合曲线也更加接近数据点。因此，在使用最小二乘进行数据拟合时应尽量使用较高次数的拟合函数。通过对比图像，可以发现二次函数和指数型二次函数图像比较接近，那是因为指数型二次函数是对数据点求对数进行拟合的，因此二者比较接近。给定的数据点在中间的位置时波动较大，造成了拟合曲线的误差也较大，但并没有影响到整体的规律。在进行最小二次编程中，使用正交化方法中间变量较多容易造成混乱，计算量也较大，但是比使用法方程方法更稳定。本题在编程过程中有一定的适用性，只需要对数据点进行更改就可以拟合出给定次数的多项式，为以后的数据拟合提供了方便。
4.设计算法，求出非线性方程52645200x x -+=的所有根，并使误差不超过410-。
问题分析：
本题采用牛顿迭代法计算方程根。可令5
2
()64520f x x x =-+ ，其迭代格式为
1'()
(1)()
k k k k f x x k x f x ?+=+=-
要使牛顿迭代法收敛，则必须寻找根的合理区间[a,b],使得在该区间内()()0f a f b <,即有根。在选定的区间内函数()f x 的一二阶导数'
()f x 、''
()f x 均不变号。选定一个初值0x ，使
''00()()0f x f x > 则牛顿迭代法收敛。
程序实现：
clear;clc;
%%%观察根的大致位置 a=-10:0.1:10; for i=1:201
y(i)=6*a(i)^5-45*a(i)^2+20; end plot(a,y); grid; k=1; y(202)=0;
%%%使用二分法来确定迭代初值 for i=1:201 if y(i)*y(i+1)<0 b(k)=a(i); k=k+1; continue end end
x=b;
n=length(x);
%%%牛顿迭代法计算方程根
dlt=1e-4;
for i=1:n
for j=1:100
X(i,j)=x(i)-(6*x(i)^5-45*x(i)^2+20)/(30*x(i)^4-90*x(i));
c=abs(X(i,j)-x(i));
if c break end x(i)=X(i,j); end end disp('计算结果为:'); x 运行结果： 函数图像如下，可以观察根的大致位置。 计算结果： 总结心得： 在该题的计算中，首先画出了函数的图像可以初步观察出根的大致位置，有效缩减了求根区间。为防止迭代进入死循环，还设置了最高迭代次数为100次，结果显示方程的三个根均在4次以内就可以迭代出，可以看出牛顿迭代法的求解速度相当快。 5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。 问题分析： 高斯消去法在整个消去过程中，(1) k kk a - 均不为零，但是若在某一步中，对角线上元的绝对值比其下方的某些元素(1) ()k ik a i k ->的绝对值小得多，即 (1)(1) ,k k kk ik a a i k --<<> 那么，相应的数(1)(1) /k k ik ik kk l a a --= 将会是绝对值很大的数，这可能引起上溢而中断计算，即使不发生上溢，也会带来很大的舍入误差。而列主元消去法的原则使(1)()k ik a i k ->中，绝 对值最大的一个换到（k,k ）位置，这样就有效的避免了大数除以小数而造成的上溢现象。 列主元消去法的主要核心是查找最大的(1)()k ik a i k ->并交换到(1)k kk a -。 程序实现： （1）非压缩格式dat141.dat;dat142.dat clc;clear; %%%读取数据 fid=fopen('dat141.dat','r'); a=fread(fid,5,'long',0); n=a(3); a=fread(fid,[n,n+1],'float'); for i=1:n for j=1:n A(i,j)=a(j,i); end b(i)=a(i,n+1); end %%%查找主元并交换位置 for k=1:n for m=k:n if abs(A(m,k))~=max(abs(A(:,k))) continue else r=A(k,:);A(k,:)=A(m,:);A(m,:)=r; r=b(k);b(k)=b(m);b(m)=r; continue end end for i=k+1:n m=A(i,k)/A(k,k); for j=1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end end %%%回代过程 x(n)=b(n)/A(n,n); b(n)=x(n); for k=n-1:(-1):1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); b(k)=x(k); end x=x' （2）压缩格式dat143.dat clear;clc; format short %%%读取数据 fid=fopen('dat143.dat','r'); d=fread(fid,5,'long',0); n=d(3);%矩阵维数 p=d(5);%上带宽 q=d(4);%下带宽 a=fread(fid,[p+q+1,n],'float'); b=fread(fid,[n,1],'float'); a=a'; d=b'; %%%消去过程 for k=1:n-1 for i=k+1:min(k+p,n) a(i,k-i+p+1)=a(i,k-i+p+1)/a(k,p+1); for j=k+1:min(k+q,n) a(i,j-i+p+1)=a(i,j-i+p+1)-a(i,k-i+p+1)*a(k,j-k+p+1); end b(i) = b(i) - b(k) * a(i,k-i+p+1); end end %%%回代过程 x(n) = b(n)/a(n,p+1); for k=n-1:-1:1 sum=0; for j=k+1:min(k+q,n) sum=sum+a(k,j-k+p+1)*x(j); end x(k)=(b(k)-sum)/a(k,p+1); end x 运算结果： 1）非压缩格式dat141.dat;dat142.dat的运算结果如下 （2）压缩格式dat143.dat 总结心得： 所给数据的方程的根都相等，其中dat141.dat ，为10阶方程，根均为1.0000；dat142.dat 为2050阶，其根均为1.9600；dat143.dat 为43512阶，其根均为3.1413。在求解非压缩格式的方程根时，采用一般列主元高斯消去法，计算量大，耗时长。在求解压缩格式的方程根时，由于阶数较大，运算时间太长，对程序进行了改进，只对非零元素进行计算，运算步数大量减少，速度得到了极大的提高。 本专业算例：一底边绝热的正方形导热物体，内无热源，其余三边温度如下，求解该正方形内节点1、2、3、4的温度。 解：导热物体无内热源，物性为常数的控制方程可以写成 02222=??+??y T x T 考虑到空间步长相同时，离散方程可以写成： 0=-+++P P S S N N W W E E T a T a T a T a T a 式中4,1,1,1,1=====P S N W E a a a a a 从而有： 04=-+++P S N W E T T T T T 所有的节点方程可整理为 704132-=-+T T T 504241-=-+T T T 153314-=-+T T T 103423-=-+T T T clear;clc; A=[-4 1 1 0 -70;1 -4 0 1 -50;1 0 -3 1 -15;0 1 1 -3 -10]; [m,n]=size(A); for i=1:(m-1) max=i; for k=i+1:m if abs(A(k,i))>abs(A(max,i)) max=k; end end temp=A(i,i:m);A(i,i:m)=A(max,i:m);A(max,i:m)=temp; for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); end end disp('回代求解') x(m)=A(m,n)/A(m,m); for k=(m-1):-1:1 x(k)=(A(k,n)-A(k,k+1:m)*x(k+1:m)')/A(k,k); end x 运算结果为T1=28.7368，T2=24.2632，T3=20.6842，T4=18.3158。运算结果和正确结果吻合。 《计算方法》实验报告 姓名： 班级： 学号： 实验日期: 2011年10月26日 一、实验题目： 数值积分 二、实验目的： 1．熟悉matlab 编写及运行数值计算程序的方法。 2．进一步理解数值积分的基础理论。 3．进一步掌握应用不同的数值积分方法求解给定的积分并给出数据结果及误差分析。 三、实验内容： 1．分别用复合梯形求积公式及复合辛普森求积公式计算积分xdx x ln 10 ? , 要求计算精度达到410-,给出计算结果并比较两种方法的计算节点数. 2．用龙贝格求积方法计算积分dx x x ?+3 021,使误差不超过510-. 3．用3=n 的高斯-勒让德公式计算积分?3 1 sin x e x ,给出计算结果. 4．用辛普森公式(取2==M N ) 计算二重积分.5 .00 5 .00 dydx e x y ? ? - 四、实验结果： 1.（1）复合梯形法: 将区间[a,b]划分为n 等份，分点n k n a b h kh a x k ,2,1,0,,=-=+=在每个区间[1,+k k x x ](k=0,1,2,···n-1)上采用梯形公式，则得 )()]()([2)()(1 11 1 f R x f x f h dx x f dx x f I n n k k k b a n k x x k k ++===∑?∑? -=+-=+ 故)]()(2)([21 1 b f x f a f h T n k k n ++=∑-=称为复合梯形公式 计算步长和划分的区间 Eps=1E-4 h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) h1 =0.0600 N1=ceil(1/h1) N1 =17 用复合梯形需要计算17个结点。 复合梯形： function T=trap(f,a,b,n) h=(b-a)/n; . / 《计算方法》上机实验报告 班级：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 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。 《计算方法》上机实验报告 班级：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 四、运行结果： 实验二： 本科实验报告 课程名称：计算机数值方法 实验项目：方程求根、线性方程组的直接解 法、线性方程组的迭代解法、代数插值和最 小二乘拟合多项式 实验地点：行勉楼 专业班级： ******** 学号： ********* 学生姓名： ******** 指导教师：李誌，崔冬华 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" 心得体会 使用不同的方法，可以不同程度的求得方程的解，通过二分法计算的程序实现更加了解二分法的特点，二分法过程简单，程序容易实现，但该方法收敛比较慢一般用于求根的初始近似值，不同的方法速度不同。面对一个复杂的问题，要学会简化处理步骤，分步骤一点一点的循序处理，只有这样，才能高效的解决一个复杂问题。 计算方法上机报告 姓名： 学号： 班级：能动上课班级： 题目及求解： 一、对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=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)体系设计： 程序是一个简单的计算器，能正确输入数据，能实现加、减、乘、除等算术运算，运算结果能正确显示，可以清楚数据等。 (2)设计思路： 1)先在Visual C++ 6.0中建立一个MFC工程文件，名为 calculator. 2)在对话框中添加适当的编辑框、按钮、静态文件、复选框和单 选框 3)设计按钮，并修改其相应的ID与Caption. 4)选择和设置各控件的单击鼠标事件。 5)为编辑框添加double类型的关联变量m_edit1. 6)在calculatorDlg.h中添加math.h头文件，然后添加public成 员。 7)打开calculatorDlg.cpp文件，在构造函数中，进行成员初始 化和完善各控件的响应函数代码。 (3)程序清单： ●添加的public成员： double tempvalue; //存储中间变量 double result; //存储显示结果的值 int sort; //判断后面是何种运算：1.加法2.减法3. 乘法 4.除法 int append; //判断后面是否添加数字 ●成员初始化： CCalculatorDlg::CCalculatorDlg(CWnd* pParent /*=NULL*/) : CDialog(CCalculatorDlg::IDD, pParent) { //{{AFX_DATA_INIT(CCalculatorDlg) m_edit1 = 0.0; //}}AFX_DATA_INIT // Note that LoadIcon does not require a subsequent DestroyIcon in Win32 m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME); tempvalue=0; result=0; sort=0; append=0; } 数值计算方法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（说明理由）． 实验报告名称 班级：学号：姓名：成绩： 1实验目的 1）通过对二分法与牛顿迭代法作编程练习与上级运算，进一步体会二分法与牛顿迭代法的不同特点。 2）编写割线迭代法的程序，求非线性迭代法的解，并与牛顿迭代法。 2 实验内容 用牛顿法和割线法求下列方程的根 x^2-e^x=0; x*e^x-1=0; lgx+x-2=0; 3实验步骤 1)根据二分法和牛顿迭代法，割线法的算法编写相应的求根函数； 2)将题中所给参数带入二分法函数，确定大致区间； 3)用牛顿迭代法和割线法分别对方程进行求解； 3 程序设计 牛顿迭代法x0=1.0; N=100; k=0; eps=5e-6; delta=1e-6; while(1) x1=x0-fc1(x0)/fc2(x0); k=k+1; if k>N disp('Newmethod failed') break end if(abs(x1-x0) 计算方法实验报告格式 小组名称: 组长姓名(班号)： 小组成员姓名(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一个完整的实验，应包括数据准备、理论基础、实验内容及方法，最终对实验结果进行分析，以达到对理论知识的感性认识，进一步加深对相关算法的理解，数值实验以实验报告形式完成，实验报告格式如下： 一、实验名称 实验者可根据报告形式需要适当写出． 二、实验目的及要求 首先要求做实验者明确，为什么要做某个实验，实验目的是什么，做完该实验应达到什么结果，在实验过程中的注意事项，实验方法对结果的影响也可以以实验目的的形式列出． 三、算法描述（实验原理与基础理论） 数值实验本身就是为了加深对基础理论及方法的理解而设置的，所以要求将实验涉及到的理论基础，算法原理详尽列出． 四、实验内容 实验内容主要包括实验的实施方案、步骤、实验数据准备、实验的算法以及可能用到的仪器设备． 五、程序流程图 画出程序实现过程的流程图，以便更好的对程序执行的过程有清楚的认识，在程序调试过程中更容易发现问题． 六、实验结果 实验结果应包括实验的原始数据、中间结果及实验的最终结果，复杂的结果可以用表格 形式列出，较为简单的结果可以与实验结果分析合并出现． 七、实验结果分析 实验结果分析包括对对算法的理解与分析、改进与建议． 数值实验报告范例 为了更好地做好数值实验并写出规范的数值实验报告，下面给出一简单范例供读者参考． 数值实验报告 小组名称: 小组成员(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一、实验名称 误差传播与算法稳定性． 二、实验目的 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，则 《数值计算方法》上机实验报告华北电力大学 实验名称数值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.用Newton 法求方程 X 7-X 4+14=0 在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据： 设函数在有限区间[a ，b]上二阶导数存在，且满足条件 {}α?上的惟一解在区间平方收敛于方程所生的迭代序列 迭代过程由则对任意初始近似值达到的一个中使是其中上不变号 在区间],[0)(3,2,1,0,) (') ()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20 )()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f a b c f x f b a x f b f x f k k k k k k ==- ==∈≤-≠>+ 令 )9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3 2 2 5 333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f 故以1.9为起点 ?? ?? ? ='- =+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。当前后两个的差<=ε时，就认为求出了近似的根。本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。 1.2 C语言程序原代码： #include 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和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 南京信息工程大学实验（实习）报告 一、实验目的: 用最小二乘法将给定的十个点拟合成三次多项式。 二、实验步骤： 用matlab编制以函数为基的多项式最小二乘拟合程序，并用于对下列数据作三次多项式最小二乘拟合(取权函数wi=1) x -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 y -2.30 -1 -0.14 -0.25 0.61 1.03 1.75 2.75 4.42 6.94 给定直线方程为：y=1/4*x3+1/2*x2+x+1 三、实验结论： 最小二乘法：通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。 一般地。当测量数据的散布图无明显的规律时，习惯上取n次代数多项式。 程序运行结果为： a = 0.9731 1.1023 0.4862 0.2238 即拟合的三次方程为：y=0.9731+1.1023x+0.4862*x2+0.2238*x3 -2.5 -2-1.5-1-0.5 00.51 1.52 2.5 -4-20246 81012 x 轴 y 轴 拟合图 离散点 y=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3 结论： 一般情况下，拟合函数使得所有的残差为零是不可能的。由图形可以看出最小二乘解决了残差的正负相互抵消的问题，使得拟合函数更加密合实验数据。 优点：曲线拟合是使拟合函数和一系列的离散点与观测值的偏差平方和达到最小。 缺点：由于计算方法简单，若要保证数据的精确度，需要大量的数据代入计算。 数值计算方法上机实习题 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 ，并用牛顿迭代法； 计算方法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． 设?+=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; 实验一——插值方法 实验学时：4 实验类型：设计 实验要求：必修 一 实验目的 通过本次上机实习，能够进一步加深对各种插值算法的理解；学会使用用三种类型的插值函数的数学模型、基本算法，结合相应软件（如VC/VB/Delphi/Matlab/JAVA/Turbo C ）编程实现数值方法的求解。并用该软件的绘图功能来显示插值函数，使其计算结果更加直观和形象化。 二 实验内容 通过程序求出插值函数的表达式是比较麻烦的，常用的方法是描出插值曲线上尽量密集的有限个采样点，并用这有限个采样点的连线，即折线，近似插值曲线。取点越密集，所得折线就越逼近理论上的插值曲线。本实验中将所取的点的横坐标存放于动态数组[]X n 中，通过插值方法计算得到的对应纵坐标存放 于动态数组[]Y n 中。 以Visual C++.Net 2005为例。 本实验将Lagrange 插值、Newton 插值和三次样条插值实现为一个C++类CInterpolation ，并在Button 单击事件中调用该类相应函数，得出插值结果并画出图像。CInterpolation 类为 class CInterpolation { public : CInterpolation();//构造函数 CInterpolation(float *x1, float *y1, int n1);//结点横坐标、纵坐标、下标上限 ~ CInterpolation();//析构函数 ………… ………… int n, N;//结点下标上限，采样点下标上限 float *x, *y, *X;//分别存放结点横坐标、结点纵坐标、采样点横坐标 float *p_H,*p_Alpha,*p_Beta,*p_a,*p_b,*p_c,*p_d,*p_m;//样条插值用到的公有指针，分别存放 i h ，i α，i β，i a ，i b ，i c ，i d 和i m }; 其中，有参数的构造函数为 CInterpolation(float *x1, float *y1, int n1) { //动态数组x1，y1中存放结点的横、纵坐标，n1是结点下标上限（即n1+1个结点） n=n1; N=x1[n]-x1[0]; X=new float [N+1]; x=new float [n+1]; y=new float [n+1];《计算方法》课内实验报告
计算方法上机实验报告
数值分析上机作业
计算方法上机实验报告
太原理工大学数值计算方法实验报告
西安交通大学计算方法B上机报告
数值分析上机题目详解
c 计算器实验报告
数值计算方法I上机实验考试题
计算方法第二章方程求根上机报告
计算方法实验报告格式
《数值计算方法》上机实验报告
数值分析上机实验报告
(完整版)哈工大-数值分析上机实验报告
计算方法上机实习题大作业(实验报告).
计算方法实验报告 拟合
数值计算方法上机实习题
计算方法B上机报告
(完整版)数值计算方法上机实习题答案
计算方法实验报告册
相关主题
文本预览