哈工大数值分析上机实验报告.docx
《哈工大数值分析上机实验报告.docx》由会员分享，可在线阅读，更多相关《哈工大数值分析上机实验报告.docx（37页珍藏版）》请在冰点文库上搜索。
哈工大数值分析上机实验报告
实验报告一
题目：
非线性程求解
摘要：
非线性程的解析解通常很难给出，因此线性程的数值解法就尤为重要。
本实验采用两种常见的求解法二分法和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*的迭代数列{xk}，这就是Newton法的思想。
当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。
另外，若将该迭代公式改进为
其中r为要求的程的根的重数，这就是改进的Newton法，当求解已知重数的程的根时，在同种条件下其收敛速度要比Newton法快的多。
程序设计：
本实验采用Matlab的M文件编写。
其中待求解的程写成function的式，如下
functiony=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;
iff12（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）;
%%%改进常数或重根数
miu=2;
%%%初始值x0
x0=input（'inputinitialvaluex0>>'）;
k=0;%迭代次数
max=100;%最大迭代次数
R=eval（subs（f,'x0','x'））;%求解f（x0），以确定初值x0时否就是解
while（abs（R）>1e-8）
x1=x0-miu*eval（subs（f,'x0','x'））/eval（subs（df,'x0','x'））;
R=x1-x0;
x0=x1;
k=k+1;
if（eval（subs（f,'x0','x'））<1e-10）;
break
end
ifk>max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值
ss=input（'mayberesultiserror,chooseanewx0,y/n?
>>','s'）;
ifstrcmp（ss,'y'）
x0=input（'inputinitialvaluex0>>'）;
k=0;
else
break
end
end
end
k;%给出迭代次数
x=x0；%给出解
结果分析和讨论：
1.用二分法计算程
在[1，2]的根。
（
下同）
计算结果为
x=1.523；
f（x）=-3.4311e-007；
k=18；
由f（x）知结果满足要求，但迭代次数比较多，法收敛速度比较慢。
2.用二分法计算程
在[1，1.5]的根。
计算结果为
x=1.180；
f（x）=2.4815e-006；
k=17；
由f（x）知结果满足要求，但迭代次数还是比较多。
3.用Newton法求解下列程
a）
x0=0.5；
计算结果为
x=0.978；
f（x）=2.0313e-016；
k=4；
由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快。
b）
x0=1；
c）
x0=0.45,x0=0.65；
当x0=0.45时，计算结果为
x=0.983；
f（x）=-8.4584e-014；
k=4；
由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快，实际上该程确实有真解x=0.5。
当x0=0.65时，计算结果为
x=0.000；
f（x）=0；
k=9；
由f（x）知结果满足要求，实际上该程确实有真解x=0.5，但迭代次数增多，实际上当取x0〉0.68时，x≈1，就变成了程的另一个解，这说明Newton法收敛与初值很有关系，有的时候甚至可能不收敛。
4.用改进的Newton法求解，有2重根，取
x0=0.55；并与3.中的c）比较结果。
当x0=0.55时，程序死循环，无法计算，也就是说不收敛。
改
时，结果收敛为
x=0.286；
f（x）=4.1127e-007；
k=16；
显然这个结果不是很好，而且也不是收敛至程的2重根上。
当x0=0.85时，结果收敛为
x=1.489；
f（x）=2.8737e-；
k=4；
这次达到了预期的结果，这说明初值的选取很重要，直接关系到法的收敛性，实际上直接用Newton法，在给定同样的条件和精度要求下，可得其迭代次数k=15，这说明改进后的Newton法法速度确实比较快。
结论：
对于二分法，只要能够保证在给定的区间有根，使能够收敛的，当时收敛的速度和给定的区间有关，二且总体上来说速度比较慢。
Newton法，收敛速度要比二分法快，但是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结果可能不时预期需要得结果。
改进的Newton法求解重根问题时，如果初值不当，可能会不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比Newton法快得多。
实验报告二
题目：
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=[ab]
x=zeros（n,1）;
%%%函数主体
fork=1:
n-1;
%%%是否进行主元选取
ifabs（A（k,k））yzhuyuan=1;
elseyzhuyuan=0;
end
ifyzhuyuan;
%%%%选主元
t=A（k,k）;
forr=k+1:
n;
ifabs（A（r,k））>abs（t）
p=r;
elsep=k;
end
end
%%%交换元素
ifp~=k;
forq=k:
n+1;
s=A（k,q）;
A（k,q）=A（p,q）;
A（p,q）=s;
end
end
end
%%%判断系数矩阵是否奇异或病态非常重
ifabs（A（k,k））disp（‘矩阵奇异，解可能不正确’）
end
%%%%计算消元,得三角阵
forr=k+1:
n;
m=A（r,k）/A（k,k）;
forq=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）;
fork=n-1:
-1:
1;
s=0;
forr=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消去法求解，并比较结果。
记Emax为求出的解代入程后的最大误差，按要求，计算结果如下：
当
时，不选主元和选主元的计算结果如下，其中前一列为不选主元结果，后一列为选主元结果，下同。
0.3910.651
2.9722.163
2.4512.721
Emax=9.2624e-，0
此时，由于
不是很小，机器误差就不是很大，由Emax可以看出不选主元的计算结果精度还可以，因此此时可以考虑不选主元以减少计算量。
当
时，不选主元和选主元的计算结果如下
1.8770.348
1.8072.174
3.7312.609
Emax=2.4668e-005，0
此时由Emax可以看出不选主元的计算精度就不好了，误差开始增大。
当
时，不选主元和选主元的计算结果如下
1.0201.000
1.6662.000
3.1110000
Emax=0.503，0
此时由Emax可以看出，不选主元的结果应该可以说是不正确了，这是由机器误差引起的。
当
时，不选主元和选主元的计算结果如下
NaN1
NaN2
NaN3
Emax=NaN,0
不选主元时，程序报错：
Warning:
Dividebyzero.。
这是因为机器计算的最小精度为10-15，所以此时的
就认为是0，故出现了错误现象。
而选主元时则没有这种现象，而且由Emax可以看出选主元时的结果应该是精确解。
结论：
采用Gauss消去法时，如果在消元时对角线上的元素始终较大（假如大于10-5），那么本法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一步，以减少机器误差带来的影响，使法得出的结果正确。
实验报告三
题目：
Rung现象产生和克服
摘要：
由于高次多项式插值不收敛，会产生Runge现象，本实验在给出具体的实例后，采用分段线性插值和三次样条插值的法有效的克服了这一现象，而且还取的很好的插值效果。
前言：
（目的和意义）
1.深刻认识多项式插值的缺点。
2.明确插值的不收敛性怎样克服。
3.明确精度与节点和插值法的关系。
数学原理：
在给定n+1个节点和相应的函数值以后构造n次的Lagrange插值多项式，实验结果表明（见后面的图）这种多项式并不是随着次数的升高对函数的逼近越来越好，这种现象就是Rung现象。
解决Rung现象的法通常有分段线性插值、三次样条插值等法。
分段线性插值：
设在区间[a,b]上，给定n+1个插值节点
a=x0和相应的函数值y0，y1，…，yn，，求作一个插值函数
，具有如下性质：
1）
，j=0，1，…，n。
2）
在每个区间[xi,xj]上是线性连续函数。
则插值函数
称为区间[a,b]上对应n个数据点的分段线性插值函数。
三次样条插值：
给定区间[a,b]一个分划
⊿：
a=x0若函数S（x）满足下列条件：
1）S（x）在每个区间[xi,xj]上是不高于3次的多项式。
2）S（x）及其2阶导数在[a,b]上连续。
则称S（x）使关于分划⊿的三次样条函数。
程序设计：
本实验采用Matlab的M文件编写。
其中待插值的程写成function的式，如下
functiony=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;
forq=1:
n+1;
l=1;%求插值基函数
fork=1:
n+1;
ifk~=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'）%作出插值函数曲线
gridon
holdon
分段线性插值源程序
clear
n=input（'将区间分为的等份数输入：
\n'）;
s=[-1+2/n*[0:
n]];%%%给定的定点，Rf为给定的函数
m=0;
hh=0.001;
forx=-1:
hh:
1;
ff=0;
fork=1:
n+1;%%%求插值基函数
switchk
case1
ifx<=s
（2）;
l=（x-s
（2））./（s
（1）-s
（2））;
else
l=0;
end
casen+1
ifx>s（n）;
l=（x-s（n））./（s（n+1）-s（n））;
else
l=0;
end
otherwise
ifx>=s（k-1）&x<=s（k）;
l=（x-s（k-1））./（s（k）-s（k-1））;
elseifx>=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'）;
gridon
holdon
三次样条插值源程序：
（采用第一边界条件）
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;
fork=1:
n;%取出节点间距
h（k）=s（k+1）-s（k）;
end
fork=1:
n-1;%求出系数向量lamuda,miu
la（k）=h（k+1）/（h（k+1）+h（k））;
miu（k）=1-la（k）;
end
%%%%赋值系数矩阵A
fork=1:
n-1;
forp=1:
n-1;
switchp
casek
A（k,p）=2;
casek-1
A（k,p）=miu（p+1）;
casek+1
A（k,p）=la（p-1）;
otherwise
A（k,p）=0;
end
end
end
%%%%求出d阵
fork=1:
n-1;
switchk
case1
d（k）=6*f2c（[s（k）s（k+1）s（k+2）]）-miu（k）*v;
casen-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;
forx=a:
hh:
b;
ifx==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'）
holdon
gridon
结果分析和讨论：
本实验采用函数
进行数值插值，插值区间为[-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的线性分段插值函数。
xf（x）L10（x）S10（x）S40（x）X10（x）X40（x）
-1.0000.1540.1540.1540.1540.1540.154
-0.0000.2391.9200.0400.2390.9100.239
-0.0000.9411.9260.4580.9410.6650.941
-0.0000.3440.9820.9790.3440.4210.344
-0.0000.1760.1760.1760.1760.1760.176
-0.0000.378-0.6740.7440.3780.8820.378
-0.0000.321-0.2500.8660.3210.5880.321
-0.0000.649-0.4180.8490.6490.2940.649
-0.0000.0000.0000.0000.0000.0000.000
-0.0000.7880.2570.7130.7880.0000.788
-0.0000.2760.1030.7300.2760.0000.276
-0.0000.8250.2670.8830.8250.0000.825
-0.0000.0000.0000.0000.0000.0000.000
-0.0000.3850.3760.4640.3850.0000.385
-0.0000.2310.0800.8600.2310.0000.231
-0.0000.9020.7890.3270.9020.0000.902
-0.0000.0000.0000.0000.0000.0000.000
-0.0000.0000.3400.4310.0000.0000.000
-0.0000.0000.8900.8280.0000.0000.000
-0.0000.8240.0730.8100.8240.0000.824
01.0001.0001.0001.0001.0001.000
0.0000.8240.0730.8100.8240.0000.824
0.0000.0000.8900.8280.0000.0000.000
0.0000.0000.3400.4310.0000.0000.000
0.0000.0000.0000.0000.0000.0000.000
0.0000.9020.7890.3270.9020.0000.902
0.0000.2310.0800.8600.2310.0000.231
0.0000.3850.3760.4640.3850.0000.385
0.0000.0000.0000.0000.0000.0000.000
0.0000.8250.2670.8830.8250.0000.825
0.0000.2760.1030.7300.2760.0000.276
0.0000.7880.2570.7130.7880.0000.788
0.0000.0000.0000.0000.0000.0000.000
0.0000.649-0.4180.8490.6490.2940.649
0.0000.321-0.2500.8660.3210.5880.321
0.0000.378-0.6740.7440.3780.8820.378
0.0000.1760.1760.1760.1760.1760.176
0.0000.3440.9820.9790.3440.4210.344
0.0000.9411.9260.4580.9410.6650.941
0.0000.2391.9200.0400.2390.9100.239
1.0000.1540.1540.1540.1540.1540.154
从以上结果可以看到，用三次样条插值和线性分段插值，不会出现多项式插值是出现的Runge现象，插值效果明显提高。
进一步说，为了提高插值精度，用三次样条插值和线性分段插值是可以增加插值节点的办法来满足要求，而用多项式插值函数时，节点数的增加必然会使多项式的次数增加，这样会引起数值不稳定，所以说这两种插值要比多项式插值好的多。
而且在给定节点数的条件下，三次样条插值的精度要优于线性分段插值，曲线的光滑性也要好一些。
实验报告四
题目：
多项式最小二乘法
摘要：
对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。
又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。
最小二乘法是求解这一问题的很好的法，本实验运用这一法实现对给定数据的拟合。
前言：
（目的和意义）
1.学习使用最小二成法的原理
2.了解法程的特性
数学原理：
对于给定的测量数据（
升级会员 