哈工大数值分析上机实验报告 20xx年
实验报告一
题目： Gauss列主元消去法
摘要：求解线性方程组的方法很多，主要分为直接法和间接法。本实验运用直接法的Guass消去法,并采用选主元的方法对方程组进行求解。
前言：（目的和意义）
1. 学习Gauss消去法的原理。
2. 了解列主元的意义。
3. 确定什么时候系数阵要选主元
数学原理：
由于一般线性方程在使用Gauss消去法求解时，从求解的过程中可以看到，若
=0，则必须进行行交换，才能使消去过程进行下去。有的时候即使
0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行r，使得
并将第r行和第k行的元素进行交换，以使得当前的的数值比0要大的多。这种列主元的消去法的主要步骤如下：
1. 消元过程
对k=1,2,…,n-1,进行如下步骤。
1) 选主元，记
若
很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。
2) 交换增广阵A的r，k两行的元素。
(j=k,…,n+1)
3) 计算消元
(i=k+1,…,n; j=k+1,……,n+1)
2. 回代过程
对k= n, n-1,…,1,进行如下计算
至此，完成了整个方程组的求解。
程序设计：
本实验采用Matlab的M文件编写。
Gauss消去法源程序：
clear
a=input('输入系数阵：>>\n')
b=input('输入列阵b：>>\n')
n=length(b);
A=[a b]
x=zeros(n,1);
%%%函数主体
for k=1:n-1;
%%%是否进行主元选取
if abs(A(k,k))<yipusilong;%事先给定的认为有必要选主元的小数
yzhuyuan=1;
else yzhuyuan=0;
end
if yzhuyuan;
%%%%选主元
t=A(k,k);
for r=k+1:n;
if abs(A(r,k))>abs(t)
p=r;
else p=k;
end
end
%%%交换元素
if p~=k;
for q=k:n+1;
s=A(k,q);
A(k,q)=A(p,q);
A(p,q)=s;
end
end
end
%%%判断系数矩阵是否奇异或病态非常严重
if abs(A(k,k))< yipusilong
disp(‘矩阵奇异，解可能不正确’)
end
%%%%计算消元,得三角阵
for r=k+1:n;
m=A(r,k)/A(k,k);
for q=k:n+1;
A(r,q)=A(r,q)-A(k,q)*m;
end
end
end
%%%%求解x
x(n)=A(n,n+1)/A(n,n);
for k=n-1:-1:1;
s=0;
for r=k+1:n;
s=s+A(k,r)*x(r);
end
t=(A(k,n+1)-s)
x(k)=(A(k,n+1)-s)/A(k,k)
end
结果分析和讨论：
例：求解方程
。
求解的结果为：
=
例：求解方程
求得的结果为：=
结论：
采用Gauss消去法时，如果在消元时对角线上的元素始终较大（假如大于10-5），那么本方法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一步，以减少机器误差带来的影响，使方法得出的结果正确。
实验报告二
题目： Rung现象产生和克服
摘要：由于高次多项式插值不收敛，会产生Runge现象，本实验在给出具体的实例后，采用分段线性插值和三次样条插值的方法有效的克服了这一现象，而且还取的很好的插值效果。
前言：（目的和意义）
1. 深刻认识多项式插值的缺点。
2. 明确插值的不收敛性怎样克服。
3. 明确精度与节点和插值方法的关系。
数学原理：
在给定n+1个节点和相应的函数值以后构造n次的Lagrange插值多项式，实验结果表明（见后面的图）这种多项式并不是随着次数的升高对函数的逼近越来越好，这种现象就是Rung现象。
解决Rung现象的方法通常有分段线性插值、三次样条插值等方法。
分段线性插值：
设在区间[a, b]上，给定n+1个插值节点
a=x0<x1<…<xn=b
和相应的函数值y0，y1，…，yn，，求作一个插值函数
，具有如下性质：
1) 
，j=0，1，…，n。
2) 在每个区间[xi, xj]上是线性连续函数。则插值函数称为区间[a, b]上对应n个数据点的分段线性插值函数。
三次样条插值：
给定区间[a, b]一个分划
⊿：a=x0<x1<…<xN=b
若函数S(x)满足下列条件：
1) S(x)在每个区间[xi, xj]上是不高于3次的多项式。
2) S(x)及其2阶导数在[a, b]上连续。则称S(x)使关于分划⊿的三次样条函数。
程序设计流程：
本实验采用Matlab的M文件编写。其中待插值的方程写成function的方式，如下
function y=f(x);
y=1/（1+25*x*x）;
写成如上形式即可，下面给出主程序
Lagrange插值源程序：
n=input('将区间分为的等份数输入：\n');
s=[-1+2/n*[0:n]];%%%给定的定点，Rf为给定的函数
x=-1:0.01:1;
f=0;
for q=1:n+1;
l=1;%求插值基函数
for k=1:n+1;
if k~=q;
l=l.*(x-s(k))./(s(q)-s(k));
else
l=l;
end
end
f=f+Rf(s(q))*l;%求插值函数
end
plot(x,f,'r')%作出插值函数曲线
grid on
hold on
分段线性插值源程序
clear
n=input('将区间分为的等份数输入：\n');
s=[-1+2/n*[0:n]];%%%给定的定点，Rf为给定的函数
m=0;
hh=0.001;
for x=-1:hh:1;
ff=0;
for k=1:n+1;%%%求插值基函数
switch k
case 1
if x<=s(2);
l=(x-s(2))./(s(1)-s(2));
else
l=0;
end
case n+1
if x>s(n);
l=(x-s(n))./(s(n+1)-s(n));
else
l=0;
end
otherwise
if x>=s(k-1)&x<=s(k);
l=(x-s(k-1))./(s(k)-s(k-1));
else if x>=s(k)&x<=s(k+1);
l=(x-s(k+1))./(s(k)-s(k+1));
else
l=0;
end
end
end
ff=ff+Rf(s(k))*l;%%求插值函数值
end
m=m+1;
f(m)=ff;
end
%%%作出曲线
x=-1:hh:1;
plot(x,f,'r');
grid on
hold on
三次样条插值源程序：（采用第一边界条件）
clear
n=input('将区间分为的等份数输入：\n');
%%%插值区间
a=-1;
b=1;
hh=0.001;%画图的步长
s=[a+(b-a)/n*[0:n]];%%%给定的定点，Rf为给定的函数
%%%%第一边界条件Rf"(-1),Rf"(1)
v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4;
for k=1:n;%取出节点间距
h(k)=s(k+1)-s(k);
end
for k=1:n-1;%求出系数向量lamuda,miu
la(k)=h(k+1)/(h(k+1)+h(k));
miu(k)=1-la(k);
end
%%%%赋值系数矩阵A
for k=1:n-1;
for p=1:n-1;
switch p
case k
A(k,p)=2;
case k-1
A(k,p)=miu(p+1);
case k+1
A(k,p)=la(p-1);
otherwise
A(k,p)=0;
end
end
end
%%%%求出d阵
for k=1:n-1;
switch k
case 1
d(k)=6*f2c([s(k) s(k+1) s(k+2)])-miu(k)*v;
case n-1
d(k)=6*f2c([s(k) s(k+1) s(k+2)])-la(k)*v;
otherwise
d(k)=6*f2c([s(k) s(k+1) s(k+2)]);
end
end
%%%%求解M阵
M=A\d';
M=[v;M;v];
%%%%
m=0;
f=0;
for x=a:hh:b;
if x==a;
p=1;
else
p=ceil((x-s(1))/((b-a)/n));
end
ff1=0;
ff2=0;
ff3=0;
ff4=0;
m=m+1;
ff1=1/h(p)*(s(p+1)-x)^3*M(p)/6;
ff2=1/h(p)*(x-s(p))^3*M(p+1)/6;
ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p));
ff4=Rf(s(p))-M(p)*h(p)*h(p)/6;
f(m)=ff1+ff2+ff3+ff4 ;
end
%%%作出插值图形
x=a:hh:b;
plot(x,f,'k')
hold on
grid on
结果分析和讨论：
本实验采用函数
进行数值插值，插值区间为[-1，1]，给定节点为
xj=-1+jh,h=0.1，j=0,…，n。下面分别给出Lagrange插值，三次样条插值，线性插值的函数曲线和数据表。图中只标出Lagrange插值的十次多项式的曲线，其它曲线没有标出，从数据表中可以看出具体的误差。
表中，L10(x)为Lagrange插值的10次多项式，S10(x)，S40(x)分别代表n=10，40的三次样条插值函数，X10(x)，X40(x)分别代表n=10，40的线性分段插值函数。
x f(x) L10(x) S10(x) S40(x) X10(x) X40(x)
-1.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154
-0.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.04244031830239
-0.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.04705882352941
-0.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.05245901639344
-0.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176
-0.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.06639004149378
-0.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.07547169811321
-0.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.08648648648649
-0.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000
-0.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.11678832116788
-0.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.13793103448276
-0.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.16494845360825
-0.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000
-0.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.24615384615385
-0.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.35000000000000 0.30769230769231
-0.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.39024390243902
-0.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000
-0.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.64000000000000
-0.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.80000000000000
-0.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.94117647058824
0 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000
0.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.94117647058824
0.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.80000000000000
0.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.64000000000000
0.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000
0.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.39024390243902
0.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.35000000000000 0.30769230769231
0.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.24615384615385
0.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000
0.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.16494845360825
0.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.13793103448276
0.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.11678832116788
0.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000
0.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.08648648648649
0.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.07547169811321
0.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.06639004149378
0.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176
0.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.05245901639344
0.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.04705882352941
0.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.04244031830239
1.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154
从以上结果可以看到，用三次样条插值和线性分段插值，不会出现多项式插值是出现的Runge现象，插值效果明显提高。进一步说，为了提高插值精度，用三次样条插值和线性分段插值是可以增加插值节点的办法来满足要求，而用多项式插值函数时，节点数的增加必然会使多项式的次数增加，这样会引起数值不稳定，所以说这两种插值要比多项式插值好的多。而且在给定节点数的条件下，三次样条插值的精度要优于线性分段插值，曲线的光滑性也要好一些。
实验报告三
题目： 多项式最小二乘法
摘要：对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。
前言：（目的和意义）
1. 学习使用最小二成法的原理
2. 了解法方程的特性
数学原理：
对于给定的测量数据(xi,fi)(i=1,2,…，n)，设函数分布为
特别的，取
为多项式
(j=0, 1,…，m)
则根据最小二乘法原理，可以构造泛函
令
(k=0, 1,…，m)
则可以得到法方程
求该解方程组，则可以得到解
，因此可得到数据的最小二乘解
程序设计：本实验采用Matlab的M文件编写。
其中多项式函数
写成function的方式，如下
function y=fai(x,j)
y=1;
for i=1:j
y=x.*y;
end
写成如上形式即可，下面给出主程序。
多项式最小二乘法源程序
clear
%%%给定测量数据点(s,f)
s=[3 4 5 6 7 8 9];
f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05];
%%%计算给定的数据点的数目
n=length(f);
%%%给定需要拟合的数据的最高次多项式的次数
m=10;
%%%程序主体
for k=0:m;
g=zeros(1,m+1);
for j=0:m;
t=0;
for i=1:n;%计算内积(fai(si),fai(si))
t=t+fai(s(i),j)*fai(s(i),k);
end
g(j+1)=t;
end
A(k+1,:)=g;%法方程的系数矩阵
t=0;
for i=1:n;%计算内积(f(si),fai(si))
t=t+f(i)*fai(s(i),k);
end
b(k+1,1)=t;
end
a=A\b%求出多项式系数
x=[s(1):0.01:s(n)]';
y=0;
for i=0:m;
y=y+a(i+1)*fai(x,i);
end
plot(x,y)%作出拟合成的多项式的曲线
grid on
hold on
plot(s,f,'rx') %在上图中标记给定的点
结果分析和讨论：
例 用最小二乘法处理下面的实验数据.
并作出
的近似分布图。
分别采用一次，二次和五次多项式来拟合数据得到相应的拟合多项式为：
y1=-0.38643+0.82750x；
y2=-1.03024+1.06893x-0.02012x2；
y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5；
分别作出它们的曲线图，图中点划线为y1曲线，实线为y2曲线，虚线为y5曲线。’x’为给定的数据点。从图中可以看出并不是多项式次数越高越好，次数高了，曲线越能给定点处和实际吻合，但别的地方就很差了。因此，本例选用一次和两次的多项式拟合应该就可以了。
实验报告四
题目： Romberg积分法
摘要：对于实际的工程积分问题，很难应用Newton-Leibnitz公式去求解。因此应用数值方法进行求解积分问题已经有着很广泛的应用，本文基于Romberg积分法来解决一类积分问题。
前言：（目的和意义）
1. 理解和掌握Romberg积分法的原理；
2. 学会使用Romberg积分法；
3. 明确Romberg积分法的收敛速度及应用时容易出现的问题。
数学原理：
考虑积分
，欲求其近似值，通常有复化的梯形公式、Simpsion公式和Cotes公式。但是给定一个精度，这些公式达到要求的速度很缓慢。如何提高收敛速度，自然是人们极为关心的课题。为此，记T1,k为将区间[a,b]进行2k等分的复化的梯形公式计算结果，记T2,k为将区间[a,b]进行2k等分的复化的Simpsion公式计算结果，记T3,k为将区间[a,b]进行2k等分的复化的Cotes公式计算结果。根据Richardson外推加速方法，可以得到收敛速度较快的Romberg积分法。其具体的计算公式为：
1. 准备初值，计算
2. 按梯形公式的递推关系，计算
3. 按Romberg积分公式计算加速值
m=2,…，k
4. 精度控制。对给定的精度
，若
则终止计算，并取
为所求结果；否则返回2重复计算，直至满足要求的精度为止。
程序设计：本实验采用Matlab的M文件编写。
其中待积分的函数写成function的方式，例如如下
function yy=f(x,y);
yy=x.^3;
写成如上形式即可，下面给出主程序
Romberg积分法源程序
%%% Romberg积分法
clear
%%%积分区间
b=3;
a=1;
%%%精度要求
R=1e-5;
%%%应用梯形公式准备初值
T(1,1)=(b-a)*(f(b)+f(a))/2;
T(1,2)=T(1,1)/2+(b-a)/2*f((b+a)/2);
T(2,1)=(4*T(1,2)-T(1,1))/(4-1);
j=2;
m=2;
%%%主程序体%%%
while(abs(T(m,1)-T(m-1,1))>R);%%%精度控制
j=j+1;
s=0;
for p=1:2^(j-2);
s=s+f(a+(2*p-1)*h/(2^(j-1)));
end
T(1,j)=T(1,j-1)/2+h*s/(2^(j-1)); %%%梯形公式应用
for m=2:j;
k=(j-m+1);
T(m,k)=((4^(m-1))*T(m-1,k+1)-T(m-1,k))/(4^(m-1)-1);
end
end
%%%给出 Romberg积分法的函数表
I=T(m,1)
结果分析和讨论：
1. 求积分
。精确解I= 24999676。
运行程序得Romberg积分法的函数表为
1.0e+007 *
4.70101520000000 3.05022950000000 2.63753307500000
2.49996760000000 2.49996760000000 0
2.49996760000000 0 0
由函数表知Romberg积分给出的结果为2.4999676*10^7，与精确没有误差，精度很高。
2. 求积分
。
直接按前面方法进行积分，会发现系统报错，出现了0为除数的现象。出现这种情况的原因就是当x=0时，被积函数分母出现0，如果用一个适当的小数
（最好不要小于程序给定的最小误差值，但不能小于机器的最大精度）来代替可以避免这个问题。本实验取
，可得函数表为：
0.92073548319659 0.93979327500190 0.94451351171417 0.94569085359489 0.94598501993743
0.94614587227034 0.94608692395160 0.94608330088846 0.94608307538495 0
0.94608299406368 0.94608305935092 0.94608306035138 0 0
0.94608306038722 0.94608306036726 0 0 0
0.94608306036718 0 0 0 0
故该函数的积分为0.94608306036718，取8位有效数字。
3. 求积分
本题的解析解很难给出，但运用Romberg积分可以很容易给出近似解，函数表为：
0.42073549240395 0.33406972582924 0.31597536075922 0.31168023948094 0.31062036680949 0.31035626065456
0.30518113697100 0.30994390573588 0.31024853238818 0.31026707591900 0.31026822526959 0
0.31026142365354 0.31026884083167 0.31026831215439 0.31026830189296 0 0
0.31026895856465 0.31026830376269 0.31026830173008 0 0 0
0.31026830119484 0.31026830172211 0 0 0 0
0.31026830172262 0 0 0 0 0
故该函数的积分为0.31026830172262，取8位有效数字。
结论：
Romberg积分通常要求被积函数在积分区间上没有奇点。如有奇点，且奇点为第一间断点，那么采用例3的方法，还是能够求出来的，否则，必须采用其它的积分方法。当然，Romberg积分的收敛速度还是比较快的。
-
哈工大数电实验报告
HarbinInstituteofTechnology数字逻辑电路与系统课程名称院系班级哈尔滨工业大学20xx年11月实验二时序逻…
-
Java程序设计实验报告 哈工大
Java语言基础实验指导报告书20xx20xx年秋季学期使用班级姓名计算机科学与技术学院实验一Java编程开发环境一实验学时1学时…
-
哈工大数据库实验报告-实验一
数据库系统应用实验报告1班号学号姓名所用实验数据表给出创建表的sql语句createtable学院学院号char2NOTNULLU…
-
哈工大微波实验报告
电磁场与微波技术实验报告HarbinInstituteofTechnology电磁场与微波技术实验报告课程名称电磁场与微波技术设计…
-
哈工大倒立摆实验报告
研究生自动控制专业实验地点：A区主楼518房间平面二级倒立摆系统实验报告哈工大航天学院控制科学实验室平面二级倒立摆控制系统实验报告…
-
分析化学实验 氢氧化钠溶液浓度的标定 实验报告
实验报告姓名班级同组人项目氢氧化钠标准溶液浓度的标定课程分析化学学号一实验目的1学会称量瓶电子天平滴定管等常用滴定仪器的准备和使用…
- 数值分析实验报告1
-
西安交大编译原理语法分析器实验报告
编译原理课内实验报告学生姓名石磊专业班级计算机26学号2120xx5140所在学院电信学院提交日期20xx122一实验目的1强化对…
-
数据分析实验报告
数据分析课程实验报告学院理学院专业信息与计算科学班级姓名学号一实验题目所做实验属于哪一部分的内容例如一元线形回归及其在SPSS中的…
-
实验15 光亮镀锌及化学镀镍实验报告
光亮镀锌及化学镀镍1实验目的11学习和实践氯化钾光亮镀锌的实验室基本操作流程了解电镀的基本原理和工艺12学习并掌握化学镀镍的原理及…
-
哈工大近代光学实验报告
近代光学创新实验双曝光全息照相技术介绍院系英才学院专业机械设计制造及其自动化学生学号班号20xx年6月双曝光全息照相技术介绍摘要双…