数值分析实验题集锦.docx
《数值分析实验题集锦.docx》由会员分享，可在线阅读，更多相关《数值分析实验题集锦.docx（26页珍藏版）》请在冰点文库上搜索。
数值分析实验题集锦
实验一
题目：
多项式最小二乘法
1.数学原理：
对于给定的测量数据（xi,fi）（i=1,2,…，n），设函数分布为
特别的，取
为多项式
（j=0,1,…，m）
则根据最小二乘法原理，可以构造泛函
令
（k=0,1,…，m）
则可以得到法方程
求该解方程组，则可以得到解
，因此可得到数据的最小二乘解
2.程序设计：
本实验采用Matlab的M文件编写。
其中多项式函数
写成function的方式，如下
functiony=fai（x,j）
y=1;
fori=1:
j
y=x.*y;
end
写成如上形式即可，下面给出主程序。
多项式最小二乘法源程序
clear
%%%给定测量数据点（s,f）
s=[3456789];
f=[2.012.983.505.025.476.027.05];
%%%计算给定的数据点的数目
n=length（f）;
%%%给定需要拟合的数据的最高次多项式的次数
m=10;
%%%程序主体
fork=0:
m;
g=zeros（1,m+1）;
forj=0:
m;
t=0;
fori=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;
fori=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;
fori=0:
m;
y=y+a（i+1）*fai（x,i）;
end
plot（x,y）%作出拟合成的多项式的曲线
gridon
holdon
plot（s,f,'rx'）%在上图中标记给定的点
3.结果分析和讨论：
例用最小二乘法处理下面的实验数据.
xi
3
4
5
6
7
8
9
fi
2.01
2.98
3.50
5.02
5.47
6.02
7.05
并作出
的近似分布图。
分别采用一次，二次和五次多项式来拟合数据得到相应的拟合多项式为：
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’为给定的数据点。
从图中可以看出并不是多项式次数越高越好，次数高了，曲线越能给定点处和实际吻合，但别的地方就很差了。
因此，本例选用一次和两次的多项式拟合应该就可以了。
实验二
题目：
非线性方程求解
1.数学原理：
对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法：
二分法和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法快的多。
2.程序设计：
本实验采用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；%给出解
3.结果分析和讨论：
1.用二分法计算方程
在[1，2]内的根。
（
下同）
计算结果为
x=1.40441513061523；
f（x）=-3.797205105904311e-007；
k=18；
由f（x）知结果满足要求，但迭代次数比较多，方法收敛速度比较慢。
2.用二分法计算方程
在[1，1.5]内的根。
计算结果为
x=1.32471847534180；
f（x）=2.209494846194815e-006；
k=17；
由f（x）知结果满足要求，但迭代次数还是比较多。
3.用Newton法求解下列方程
a）
x0=0.5；
计算结果为
x=0.56714329040978；
f（x）=2.220446049250313e-016；
k=4；
由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快。
b）
x0=1；
c）
x0=0.45,x0=0.65；
当x0=0.45时，计算结果为
x=0.49999999999983；
f（x）=-8.362754932994584e-014；
k=4；
由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快，实际上该方程确实有真解x=0.5。
当x0=0.65时，计算结果为
x=0.50000000000000；
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.50000087704286；
f（x）=4.385198907621127e-007；
k=16；
显然这个结果不是很好，而且也不是收敛至方程的2重根上。
当x0=0.85时，结果收敛为
x=1.00000000000489；
f（x）=2.394337647718737e-023；
k=4；
这次达到了预期的结果，这说明初值的选取很重要，直接关系到方法的收敛性，实际上直接用Newton法，在给定同样的条件和精度要求下，可得其迭代次数k=15，这说明改进后的Newton法法速度确实比较快。
4.结论：
对于二分法，只要能够保证在给定的区间内有根，使能够收敛的，当时收敛的速度和给定的区间有关，二且总体上来说速度比较慢。
Newton法，收敛速度要比二分法快，但是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结果可能不时预期需要得结果。
改进的Newton法求解重根问题时，如果初值不当，可能会不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比Newton法快得多。
实验三
题目：
Gauss列主元消去法
1.数学原理：
由于一般线性方程在使用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,进行如下计算
至此，完成了整个方程组的求解。
2.程序设计：
本实验采用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
3.结果分析和讨论：
例：
求解方程
。
其中
为一小数，当
时，分别采用列主元和不列主元的Gauss消去法求解，并比较结果。
记Emax为求出的解代入方程后的最大误差，按要求，计算结果如下：
当
时，不选主元和选主元的计算结果如下，其中前一列为不选主元结果，后一列为选主元结果，下同。
0.999999347683910.99999934782651
2.000002174219722.00000217391163
2.999997608594512.99999760869721
Emax=9.301857062382624e-010，0
此时，由于
不是很小，机器误差就不是很大，由Emax可以看出不选主元的计算结果精度还可以，因此此时可以考虑不选主元以减少计算量。
当
时，不选主元和选主元的计算结果如下
1.000017846308770.99999999999348
1.999980097208072.00000000002174
3.000006634247312.99999999997609
Emax=2.036758973744668e-005，0
此时由Emax可以看出不选主元的计算精度就不好了，误差开始增大。
当
时，不选主元和选主元的计算结果如下
1.421085471520201.00000000000000
1.666666666666662.00000000000000
3.11111111111111300000000000000
Emax=0.70770085900503，0
此时由Emax可以看出，不选主元的结果应该可以说是不正确了，这是由机器误差引起的。
当
时，不选主元和选主元的计算结果如下
NaN1
NaN2
NaN3
Emax=NaN,0
不选主元时，程序报错：
Warning:
Dividebyzero.。
这是因为机器计算的最小精度为10-15，所以此时的
就认为是0，故出现了错误现象。
而选主元时则没有这种现象，而且由Emax可以看出选主元时的结果应该是精确解。
4.结论：
采用Gauss消去法时，如果在消元时对角线上的元素始终较大（假如大于10-5），那么本方法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一步，以减少机器误差带来的影响，使方法得出的结果正确。
实验四
1.题目
已知f（x）的观测数据
构造插值多项式。
2.算法原理
对于给定n+1个节点，n次拉格朗日插值多项式由下式给出：
其中
由于（）klx是一个关于x的n次多项式，所以L（x）为关于x的不高于n次的代数多
项式。
当ix=x时，（）iiLx=y，满足插值条件。
3.流程图
4.输出结果
x=2.5000000处的函数值为：
y=-6.3750000
5.结果分析
由于所知的拉格朗日计算机算法只能实现计算某一特定值的近似函数值，而
不知如何导出表达式，故例求x=2.5处的函数值以说明表达式以得出，只是在计
算机程序中。
并且也能达到拉格朗日插值法使用的目的。
实验题五
1.题目
用列主元消去法求解方程组
并求出系数矩阵A的行列式的值detA。
2.算法原理
2.1顺序高斯消去法
顺序高斯消去法是利用线性方程组初等变换中的一种变换，即用一个不为零
的数乘一个方程后加至另一个方程，使方程组变成同解的上三角方程组，然后再
自下而上对上三角方程组求解。
这样，顺序高斯消去法可分成“消去”和“回代”
两个过程。
在用顺序高斯消去法时，在消元之前检查方程组的系数矩阵的顺序主子式，
当阶数较高时是很难做到的。
若线性方程组的系数具有某种性质时，如常遇到的
对角占优方程组，自然能够用高斯消去法求解。
2.2列选主元消去法
线性方程组只要系数矩阵非奇异，就存在惟一解，但是按顺序消元过程中可
能出现主元素akk（k）=0，这时尽管系数矩阵非奇异，消元过程无法再进行，或者
即使akk（k）≠0，但如果其绝对值很小，用它作除数也会导致其他元素的数量级急剧增大和使舍入误差扩大，将严重影响计算的精度。
为避免在校园过程确定乘数时的所用除数是零或绝对值小的数，即零主元或小主元，在每一次消元之前，要增加一个选主元的过程，将绝对值大的元素交换到主对角线的位置上来。
列选主元是当高斯消元到第k步时，从k列的akk以下（包括akk）的各元素中选出绝对值最大的，然后通过行交换将其交换到kka的位置上。
交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。
列选主元消去法常用来求行列式。
设有矩阵
用列主元消去法将其化为上三角形矩阵，对角线上元素为
于是行列式
其中m为所进行的行交换次数。
这是实际中求行列式值的可靠方法。
3.流程图
4.输出结果
x
（1）=1.0000
X
（2）=2.0000
X（3）=3.0000
detA=-66.0000
5.结果分析
采用计算机运算在计算大数据时有明显的优点，另外也需要考虑到存储。
高斯消去法的使用条件是
而列选主元法可以保证这一
条件。
并且可以避免在消元过程确定乘数时所用除数是绝对值小的数，相对全选
主元的运算量小，一般也可以满足精度要求。
实验题六
题目：
常微分方程初值问题数值解法
一、问题提出
科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题：
（1）
分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解
。
（2）
用r=3的Adams显式和预-校式求解
取步长h=0.1，用四阶标准R-K方法求值。
（3）
用改进Euler法或四阶标准R-K方法求解
取步长0.01,计算
数值解，参考结果
。
（4）利用四阶标准R-K方法求二阶方程初值问题的数值解
（I）
（II）
（III）
（IV）
二、要求
1、根据初值问题数值算法，分别选择二个初值问题编程计算；
2、试分别取不同步长，考察某节点处
数值解的误差变化情况；
3、试用不同算法求解某初值问题，结果有何异常；
4、分析各个算法的优缺点。
三、实验程序及分析
（Ⅰ）
（1）、算法程序
functionE=Euler_1（fun,x0,y0,xN,N）
%Euler向前公式，其中
%fun为一阶微分方程的函数
%x0,y0为初始条件
%xN为取值范围的一个端点
%h为区间步长
%N为区间个数
%x为Xn构成的向量
%y为yn构成的向量
x=zeros（1,N+1）;y=zeros（1,N+1）;
x
（1）=x0;y
（1）=y0;
h=（xN-x0）/N;
forn=1:
N
x（n+1）=x（n）+h;
y（n+1）=y（n）+h*feval（fun,x（n）,y（n））;
end
T=[x',y']
functionz=f（x,y）
z=4*x/y-x*y;
（2）、运行程序
>>Euler_1（'f',0,3,2,20）
四、实验结果
>>Euler_1（'f',0,3,2,20）
T=
03.0000
0.10002.9836
0.20002.9517
0.30002.9058
0.40002.8481
0.50002.7810
0.60002.7073
0.70002.6297
0.80002.5511
0.90002.4739
1.00002.4004
1.10002.3325
1.20002.2714
1.30002.2177
1.40002.1717
1.50002.1332
1.60002.1017
1.70002.0765
1.80002.0567
1.90002.0414
2.00002.0299
五、实验结论：
欧拉（Euler）算法是数值求解中最基本、最简单的方法，但其求解精度较低。
改进欧拉格式较欧拉格式提高了精度，其截断误差比欧拉格式提高了一阶。
龙格－库塔法是用于模拟长微分方程的解的重要的一类隐式或显式迭代法。
工程中很多的地方用到龙格库塔求解微分方程的数值解，龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。
升级会员 