计算方法上机题
1. 计算100000
21
1k S k
==
∑
,要求误差小于6
10-,给出实现算法。 算法组织：利用绝对值的递增相加可以减小误差 具体算法：
%% 公式求和 %% S = 0; k = 100000;
%%利用for 循环对公式进行求和 for n = 1: 100000 S = S+1/(k*k); k = k-1; end
%%精确到小数点后6位,以保证误差范围 vpa(S,7)
运行程序后，得到的结果是1.644924，精确到了小数点后6位。
2. 编写实现对N 阶非奇矩阵A 进行LU 分解的程序。 算法组织：
（1） 根据定理可知若N 阶矩阵的各阶顺序主子式行列式不为零，则存在唯一的单位下三
角矩阵L 和上三角矩阵U ，满足A=LU ，因此首先需要判断A 中的元素a kk （k=0,1,2…，n ）不为零。
（2） 根据如下的式子计算L 和U 中的元素
121231(,,...,1,0,0...0)000j j jj ij i i i ii u u u a l l l l -?? ? ? ? ? ?= ? ? ? ? ? ???
1
1(,1,2,3
,)i ij ik kj ij k a l u u j i i n -==+≥=∑
1
(,1,2,3
,)j ij ik kj k a l u j i i n ==<=∑
当i=1时，得到u 1j =a 1j (j=1,2,3,…,n),而当j=1时，得l i1= a i1/u 11，由此可以计算出L 的第一列元素和U 的第一行元素。然后可以依据下式计算L 和U 中的其他元素。
1
1
(,1,2,3
,)i ij ij it tj t u a l u j i i i i n -==-=+++∑
1
1
1
()(1,2,3
,)i ki ki it tj t ii l a l u k i i i n u -==-=+++∑
具体算法：
[n,n]=size(A); %检测A 的阶数 L=eye(n); %初始化单位下三角阵L U=zeros(n,n); %初始化上三角阵U
for i=1:n
U(1,i)=A(1,i); %给上三角阵U 的第一行赋值
L(i,1)=A(i,1)/U(1,1); %给单位下三角阵L 的第一列赋值 end for i=2:n for j=i:n for k=1:i-1
M(k)=L(i,k)*U(k,j); end
U(i,j)=A(i,j)-sum(M); %给上三角阵U 的第i 行赋值 end
for j=i+1:n for k=1:i-1
M(k)=L(j,k)*U(k,i); end
L(j,i)=(A(j,i)-sum(M))/U(i,i); %给单位下三角阵L 的第i 列赋值 end end
%分别输出L,U 的矩阵 L U
3.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。具体的要求参见所附的相关文档。 针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。 具体算法：
fid = fopen('fun001.dat','r'); % 获取数据文件文件头 headfile = fread(fid,3,'uint32'); % 确认数据文件标识号是否正确 id = headfile(1);
if id ~= hex2dec('F1E1D1A0') fprintf('不是正确的数据文件.\n'); return ; end
% 获取数据文件版本号 ver = headfile(2); if ver == hex2dec('101') % 获取方程阶数
n = headfile(3)
% 获取稀疏矩阵
A = fread(fid,[n,n],'float32');
% 获取右端系数
B = fread(fid,[n,1],'float32');
% 对系数矩阵A和右端向量的方程组使用高斯方法求解 x=GAUSSPP(A,B,n);
end
if ver == hex2dec('102')
% 系数矩阵为条状对角阵
% 获取方程阶数
n = headfile(3)
% 获取稀疏矩阵
A = fread(fid,[n,n],'float32');
% 获取右端系数
B = fread(fid,[n,1],'float32');
% 对系数矩阵A和右端向量的方程组使用高斯方法求解x=GAUSSPP(A,B,n);
end
if ver == hex2dec('201')
% 系数矩阵为压缩格式下的条状对角阵
% 获取方程阶数
n = headfile(3)
% 获取条状矩阵的上下带宽
headfile2 = fread(fid,2,'uint32');
p = headfile2(1);
q = headfile2(2);
% 压缩系数矩阵
A = fread(fid,[p+q+1,n],'float32');
% 获取右端系数
B = fread(fid,[n,1],'float32');
% 对系数矩阵A和右端向量的方程组使用高斯方法求解x=GAUSSPP(A,B,n,p,q);
end
fclose(fid);
GAUSSPP(A,B,n)的程序如下：
function x=GAUSSPP(A,B,n)
for k=1:n-1
for i=k+1:n
l(i,k)=A(i,k)/A(k,k);
for j=k:n
A(i,j)=A(i,j)-l(i,k)*A(k,j);
end
B(i)=B(i)-l(i,k)*B(k);
end end
x(n)=B(n)/A(n,n); for k=n-1:-1:1 s=0; for j=k+1:n
s=s+A(k,j)*x(j); end
x(k)=(B(k)-s)/A(k,k); end
GAUSSPP(A,B,n,p,q)的程序如下： function x=GAUSSPP(A,B,n,p,q) for k=1:n-1
for i=k+1:min(k+p,n)
l(i,p-i+k+1)=A(i,p-i+k+1)/A(k,p+1); for j=0:q
A(i,p-i+k+1+j)=A(i,p-i+k+1+j)-l(i,p-i+k+1)*A(k,p+1+j); end
B(i)=B(i)-l(i,p-i+k+1)*B(k); end end
x(n)=B(n)/A(n,p+1); for k=n-1:-1:1 s=0;
for j=1:min(q,n-k) s=s+A(k,1+p+j)*x(j+k); end
x(k)=(B(k)-s)/A(k,p+1); end
4.已知某产品从1900年到2010年每隔10年的产量,用多项式插值和三次样条插值的方法，画出每隔一年的插值曲线的图形, 试计算并比较在不同方法下的2005年以及2015年的产量。
算法组织： 1. 多项式插值
利用多项式插值进行数据近似时要选取一个函数g k (x),通常选取如下的函数：
(),0,1,2,
,k k g x x k n ==
因此，欲寻求的函数就是n 次多项式
()k k P x a x =∑
根据差值条件可以得到下式：
2
012n
i i i n i y a a x a x a x =++++
我们可以把这个方程抽象成一个线性方程组：
2
10000
021*******n x x x a y n a y x x x a y n n m x
x x m m m ?????? ? ? ?
?
? ? ?= ? ? ?
?
? ? ? ? ? ? ????? ??
?
这样就把求解数据近似的问题转化成了求解线性方程组的问题。同时观察上式可以看出方程的系数是一个Vandermonde 矩阵。 具体算法：
% 多项式插值近似 %
Y = 1900:10:2010; A = zeros(12); YB = 1900;
% 根据数据构建11阶Vandermonde 矩阵A
for i = 1 : 12; for j = 1 : 12; A(i,j) = YB^(j-1); end
YB = YB + 10; end
% 根据数据构建矩阵B
B = [75.995;91.972;105.711;123.203;131.699;150.697;179.323;203.212; 226.505;251.525;291.854;325.433];
% 调用列主元Gauss 消去法函数,求多项式系数矩阵
x = Gauss_pivot(A,B); a =zeros(1,12); for n = 1:12 a(1,n) = x(13-n,1); end
e1 = polyval(a,[2005,2015])
T = 1900:2010
plot(Y,B,'r') % 红色表示原数据图像 hold on
plot(T,polyval(a,T),'g') % 绿色表示多项式插值近似图像 hold on
最终得到的结果是 e1 =
342 456
即2005年的产量为342,2015年的产量为456。 多项式插值的结果如图（1）所示：
图（1）多项式插值的图形
图中红色曲线表示原始数据，绿色曲线表示多项式插值的数据。
2.三次样条函数插值
算法组织：
由于分段多项式插值不够光滑，因此需要较高次的分段多项式，常用三项多项式，即在每一个子区间上市一个三次多项式，这样的分段三次多项式就是三次样条函数。
具体算法:
n = 12;
%构造数据矩阵Y
Y = 1900:10:2010;
%构造数据矩阵B
B = [75.995;91.972;105.711;123.203;131.699;150.697;179.323;203.212;
226.505;251.525;291.854;325.433];
% 构造步长矩阵h
h = zeros(1,n-1);
for k = 1:11
h(k) = Y(k+1)-Y(k);
end
% 构造M的系数矩阵A
A = zeros(n);
% 由第三类边界条件可得
A(1,2) = -2;
% 由第三类边界条件可得
A(n,n-1) = -2;
for i = 1:n
A(i,i) = 2;
end for i =2:n-1
A(i,i+1) = h(i)/(h(i-1)+h(i)); A(i,i-1) = h(i-1)/(h(i-1)+h(i)); end
% 构造M 的结果矩阵d d =zeros(n,1);
d(1) = -12*h(1)*ChaShang_4(Y,B,1,2,3,4);
d(n) = 12*h(n-1)*ChaShang_4(Y,B,n-3,n-2,n-1,n); for j = 2:n-1
d(j) = 6*ChaShang_3(Y,B,j-1,j,j+1); end
M = Gauss_pivot(A,d);
% 构造在不同区间内三次样条拟合结果的矩阵S T = 1900:2010; S = zeros(1,(n-1)*10+1); for i = 1:n-1 num = 1; for x = 1:10 S(1,(i-1)*10+num) =
M(i)*((Y(i+1)-T((i-1)*10+num))^3)/(6*h(i))+M(i+1)*((T((i-1)*10+num)-Y (i))^3)/(6*h(i))+(B(i)-((M(i)*(h(i)^2))/6))*(Y(i+1)-T((i-1)*10+num))/ h(i)+(B(i+1)-((M(i+1)*(h(i)^2))/6))*(T((i-1)*10+num)-Y(i))/h(i); num = num+1; end end
% 给出S 的边界条件
S(1,111)=M(11)*((Y(12)-T(111))^3)/(6*h(11))+M(12)*((T(111)-Y(11))^3)/(6*h(11))+(B(11)-((M(11)*(h(11)^2))/6))*(Y(12)-T(111))/h(11)+(B(12)-((M(12)*(h(11)^2))/6))*(T(111)-Y(11))/h(11); plot(T,S,'m'); hold on
5.假定某天的气温变化记录如下表所示，试用最小二乘法找出这一天的气温变化的规律,试计算这一天的平均气温。
考虑使用二次函数、三次函数、四次函数以及指数型的函数,计算相应的系
数，估算误差，并作图比较各种函数之间的区别。 算法组织：
给定数据点{（x i ,y i ）}(i=1,2,3…,m)和一组函数g k （x ），求出一系列数a 1,a 2,a 3,…,a n ,使函数1122()()()()n n p x a g x a g x a g x =++
+满足
21/22([()])i i E p x y =-∑
可以将上式，改写为如下的形式：2
2
2
21
1
[()][()]
min m n
i
i
k
k
i
i
i k E p x y a g x y ===-=-=∑∑∑
也就是要求出使上式成立的a k 的值。
因此得到了使上式成立的必要条件：2
20k
E a ?=?
最终可以得到以下的一个方程组：
1112111212222212
n n n n nn n n a b a b a b ααααααααα??????
??? ?
??? ?= ??? ?
??? ???????
解这个方程就可以得到最小二乘法的系数。 具体算法：
% 根据采样数据组构建矩阵
A2 = 0:24;
B =[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20,18 ,17,16]; x = 0:0.1:24;
% 构建最小二乘拟合法的系数矩阵A
A = zeros(25,6); for i = 1 : 25; for j = 1 : 6; A(i,j) = A2(i)^(j-1); end end
x1 = (A'*A)\(A'*B');
a1 =zeros(1,6); % 初始化系数矩阵a1
% 将系数矩阵转化为标准矩阵形式
for n = 1:6 a1(1,n) = x1(7-n,1); end
plot(x,polyval(a1,x),'b') hold on
最终的拟合结果如图(3)所示：
图（3）最小二乘法拟合结果
从图中可以得出结论，当天气温变化规律0-5 时气温略有小波动,但基本保持不变,并于4 时达到当天气温最低值; 5-15 时气温不断升高,并于15 时达到当天气温最高值; 15-24 时气温不断降低。
利用不同的函数进行拟合，拟合的结果如下：
二次函数拟合、三次函数拟合和四次函数拟合：
具体算法：
x=[0:1:24];
y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20, 18,17,16];
a1=polyfit(x,y,2) % 二次多项式拟合
a2= polyfit(x,y,3) % 三次多项式拟合
a3= polyfit(x,y,4) % 四次多项式拟合
b1= polyval(a1,x)
b2= polyval(a2,x)
b3= polyval(a3,x)
r1= sum((y-b1).^2) % 二次多项式的误差平方和
r2= sum((y-b2).^2) % 三次多项式的误差平方和
r3= sum((y-b3).^2) % 四次多项式的误差平方和
plot(x,y,'*')
hold on
plot(x,b1, 'r') % 二次多项式的图像
hold on
plot(x,b2, 'g') % 三次多项式的图像
hold on
plot(x,b3, 'b') % 四次多项式的图像
拟合图像如下图（4）所示：
图（4 ）多项式拟合的结果
r1=442.8523,r2=254.9607 ,r3 =170.5871。从误差上可以看出,多项式拟合的次数越高误差越小。
最终，拟合而成的多项式为：
二次多项式：
Linear model Poly2:
f(x) = p1*x^2 + p2*x + p3
Coefficients (with 95% confidence bounds):
p1 =-0.1 (-0.1401,-0.05989)
p2 =2.775 (1.779,3.772)
p3 =7.815 (2.651,12.98)
2
=-++
0.1 2.7757.815
y x x
三次多项式：
Linear model Poly3:
f(x) = p1*x^3 + p2*x^2 + p3*x + p4
Coefficients (with 95% confidence bounds):
p1 =-0.009389 (-0.01435, -0.004426)
p2 =0.238 (0.05662, 0.4194)
p3 =-0.4038 (-2.255, 1.447)
p4 =13.52 (8.491, 18.54)
32
=-+-+
0.0093890.2380.403813.52
y x x x
四次多项式：
Linear model Poly4:
f(x) = p1*x^4 + p2*x^3 + p3*x^2 + p4*x + p5
Coefficients (with 95% confidence bounds):
p1 =0.001012 (0.0003408, 0.001683)
p2 =-0.05796 (-0.09044, -0.02548)
p3 = 0.9777 (0.464, 1.491)
p4 = -4.168 (-7.11, -1.226)
p5 =17.2 (12.32, 22.08)
432
=-+-+
y x x x x
0.0010120.057960.9777 4.16817.2
指数函数拟合具体算法：
x=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]; y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20, 18,17,16];
opts = fitoptions('Method','Nonlinear','Normalize','On');
ftype = fittype('a*exp(-b*(x-c).^2)','options',opts);
[fresult,gof] = fit(x',y',ftype)
plot( x, fresult(x), x, y, '* ')
拟合的结果如图（5）所示：
图（5）指数拟合的结果
指数拟合的参数为：
General model:
fresult(x) = a*exp(-b*(x-c).^2)
where x is normalized by mean 12 and std 7.36
Coefficients (with 95% confidence bounds):
a=28.98 (26.34, 31.62)
b=0.3431 (0.2337, 0.4526)
c=0.3008 (0.1552, 0.4463)
所以最终的拟合而成的指数函数为：
2
(0.3431(0.3008))
=
y e-?-
28.98x
从两幅图的比较可以看出四次多项式的拟合效果比较理想。
计算方法上机实验报告
. / 《计算方法》上机实验报告 班级：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 四、运行结果： 实验二：
计算方法上机题答案
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 ，并用牛顿迭代法； 数值计算方法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． 设?+=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; 《数值计算方法》上机实验报告华北电力大学 实验名称数值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 实习题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; 实验一——插值方法 实验学时：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上机报告
数值分析上机题目详解
数值计算方法上机实习题
数值计算方法I上机实验考试题
计算方法第二章方程求根上机报告
(完整版)数值计算方法上机实习题答案
《数值计算方法》上机实验报告
数值分析上机实验报告
计算方法与实习上机题答案
太原理工大学数值计算方法实验报告
(完整版)哈工大-数值分析上机实验报告
计算方法上机实习题大作业(实验报告).
计算方法实验报告格式
计算方法上机答案
计算方法B上机报告
计算方法实验报告册