山建大数值分析实验报告Word文档下载推荐.docx
《山建大数值分析实验报告Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《山建大数值分析实验报告Word文档下载推荐.docx（30页珍藏版）》请在冰豆网上搜索。
fori=k+1:
n
a（i,k）=a（i,k）/a（k,k）;
b（i）=b（i）-a（i,k）*b（k）;
forj=k+1:
ifa（k,k）==0
主元素为零，消去法无法继续'
）;
break;
else
a（i,j）=a（i,j）-a（i,k）*a（k,j）;
end
b（n）=b（n）/a（n,n）;
fori=（n-1）:
-1:
1
w=0;
forj=（i+1）:
w=w+a（i,j）*b（j）;
b（i）=（b（i）-w）/a（i,i）;
y=b;
（b）高斯列主元消去法
functionz=gauss2（a,b,ep）
%高斯列主元素消元法
ifnargin==2
ep=0.000001
p=a（k,k）;
I=k;
fori=k:
ifabs（a（i,k））>
abs（p）
p=a（i,k）;
I=i;
ifp<
=ep
z=0;
ifI~=k
forj=k:
w=a（k,j）;
a（k,j）=a（I,j）;
a（I,j）=w;
u=b（k）;
b（k）=b（I）;
b（I）=u;
%disp（'
方程组的解为'
）;
z=b;
（5）实验结果：
>
2-100
-1-3-20
-13-20
00-35
b=
6
1
0
gauss1（a,b）
方程组的解为
ans=
35/12
-1/6
-41/24
-33/40
（6）实验体会：
主元消去法和高斯消去法的确是两个非常锻炼人编程的方法，在编写程序时，需要使用的大量的循环和分支结构，但无论是高斯消去法还是高斯列主元法，它们的原理还算不难理解，通过变成能够较好的理解它们。
实验二解线性方程组的迭代法
2.1实验目的
掌握解线性方程组的雅可比迭代和高斯-塞德尔迭代算法；
培养编程与上机调试能力.
2.2实验要求：
（1）选择一种计算机语言（Matlab）设计出雅可比（Jacobi）Gauss-Seidel、SOR迭代法，迭代法的算法程序，并且选择不同的迭代次数，观察输出结果；
（2）利用Matlab求方程组
2.3实验内容
计算书上的习题
P61例2.5.1
a=[2023;
181;
2-315];
b=[24;
12;
30];
x0=[0;
0;
0];
（2）对应程序：
Jacobi迭代法：
functionX=jacobi（a,b,X0,ep）
%Jacobi迭代法求解方程组
ifnargin==3
ep=1.0e-6;
elseifnargin<
3
error
return
D=diag（diag（a））;
D=inv（D）;
L=tril（a,-1）;
U=triu（a,1）;
B=-D*（L+U）;
f=D*b;
X=B*X0+f;
k=1;
while（norm（X-X0,inf）>
=ep）&
（k<
=1000）
X0=X;
X=B*X0+f;
k=k+1;
disp（'
迭代次数为'
k
return
高斯-塞德尔（Gauss-Seidel）：
functionz=gauss2（a,b,ep）
gauss2（a,b,ep）
（3）实验结果：
迭代次数为
k=
6
0.33
1.33
2.67
高斯消元法：
0.15
1.94
2.37
（4）实验体会：
Jacobi迭代法和高斯消元法都能很好的解决方程组的求解问题，在上机程中遇到了也遇到了不少问题，但最后在老师的悉心辅导下都得到了很好的解答，这两个程序使我明白了要变出好的程序就需要我们积极思考问题，勇于发现和解决问题。
实验三矩阵特征值问题计算
3.1实验目的
掌握求矩阵的特征值和主特征向量的幂法；
3.2实验要求
（1）选择一种计算机语言设计出幂法求主特征值和相应特征向量的程序，并且选择不同的初值，观察所需的迭代次数和迭代结果.
（2）利用Matlab求特征值和特征向量
调用格式1：
eig（A）%得到特征值列向量
调用格式2：
，其中为由特征列向量构成的方阵，为由特征值构成的对角阵.%得到特征值和所对应的特征向量
3.3实验内容
P81例3.1.1
A=[246;
3915;
41636]x0=[1;
1;
1]
乘幂法
functiony=chm（a,x0,k）
%乘幂法求主特征值及特征向量
ifnargin<
3
k=10
fori=1:
[mj]=max（abs（x0））;
y1=（1/x0（j））*x0;
%标准化
x0=a*y1;
[mj]=max（abs（x0））
lanmda=x0（j）
u1=y1
41636]
x0=[1;
chm（A,x0,6）
经典Jacobi法
functionyy=Jeig（a,eps）
%经典Jacobi法-求全部特征值及特征向量
ifa~=a'
输入错误,应输入对称矩阵'
n=size（a）;
ifnargin==1
eps=1.0e-6;
I=eye（n）;
SS=0;
fori=2:
forj=1:
i-1
SS=SS+a（i,j）^2;
end
k=1;
G=a;
R1=I;
p=0;
q=0;
while（SS>
=eps）&
10000）
a=G;
R=R1;
fori=2:
ifabs（a（i,j））>
w
w=abs（a（i,j））;
p=i;
q=j;
ifa（p,p）==a（q,q）
theta=sign（a（p,q））*pi/4;
G（p,p）=a（p,p）*（cos（theta））^2+a（q,q）*（sin（theta））^2+a（p,q）*sin（2*theta）;
G（q,q）=a（p,p）*（sin（theta））^2+a（q,q）*（cos（theta））^2-a（p,q）*sin（2*theta）;
G（p,q）=（a（q,q）-a（p,p））*（cos（theta））*（sin（theta））+a（p,q）*（cos（2*theta））;
G（q,p）=G（p,q）;
fori=1:
R1（i,p）=R（i,p）*cos（theta）+R（i,q）*sin（theta）;
R1（i,q）=-R（i,p）*sin（theta）+R（i,q）*cos（theta）;
if（i~=p）&
（i~=q）
G（p,i）=a（i,p）*（cos（theta））+a（i,q）*（sin（theta））;
G（i,p）=G（p,i）;
G（q,i）=-a（i,p）*（sin（theta））+a（i,q）*（cos（theta））;
G（i,q）=G（q,i）;
else
gasi=（a（p,p）-a（q,q））/（2*a（p,q））;
t=sign（gasi）*（-abs（gasi）+sqrt（1+gasi^2））;
G（p,p）=a（p,p）*（1/（1+t^2））+a（q,q）*（t^2/（1+t^2））+a（p,q）*2*（t/（1+t^2））;
G（q,q）=a（p,p）*（t^2/（1+t^2））+a（q,q）*（1/（1+t^2））-a（p,q）*2*（t/（1+t^2））;
G（p,q）=（a（q,q）-a（p,p））*（t/（1+t^2））+a（p,q）*（（1/（1+t^2））-（t^2/（1+t^2）））;
R1（i,p）=R（i,p）*（1/sqrt（1+t^2））+R（i,q）*（t/sqrt（1+t^2））;
R1（i,q）=-R（i,p）*（t/sqrt（1+t^2））+R（i,q）*（1/sqrt（1+t^2））;
G（p,i）=a（i,p）*（1/sqrt（1+t^2））+a（i,q）*（t/sqrt（1+t^2））;
G（q,i）=-a（i,p）*（t/sqrt（1+t^2））+a（i,q）*（1/sqrt（1+t^2））;
SS=SS+G（i,j）^2;
G的对角线即为近似的特征值'
G
R的列向量为相应的近似特征向量'
）,
R
Jeig（A,10e-9）
乘幂法：
m=
43.56
j=
3
lanmda=
u1=
0.50
0.10
1.00
经典Jacobi法：
G的对角线即为近似的特征值
G=
0.120.200.20
0.202.88-0.75
0.20-0.7544.99
R的列向量为相应的近似特征向量
R=
0.540.010.31
-0.620.210.06
0.92-0.390.44
本次程序是由老师提供的，但是同学们都认真阅读过了，我发现此程序是相当复杂的，要对矩阵进行迭代等很多操作，我还了解到我们在处理较大的问题是必须要使用到程序，因此让我对程序产生的浓厚的兴趣，同时我们也认识到我们所编写的程序稳定性很差，因此我们还需要更多的练习。
实验四插值法
4.1实验目的
1掌握插值法的基本思路和步骤；
2培养编程与上机调试能力。
4.2实验要求
用Matlab和插值中的某种具体算法编写代码并执行，完成解决具体问题。
4.3实验内容
MatlabSpline Tools
x=00.20000.40000.60000.8000
y=1.00001.22141.49181.82212.2255
x0=0.1500
functions=Newton（x,y,x0,nn）
%nn为Newton插值多项式的次数
nx=length（x）;
ny=length（y）;
ifnx~=ny
warning（'
向量x与y的长度应该相同'
return;
m=length（x0）;
m
yy=y;
t=0;
kk=1;
while（kk<
=nn）
kk=kk+1;
fork=kk:
nx
yy（k）=（yy（k）-yy（kk-1））/（x（k）-x（kk-1））;
t=yy
（1）;
fork=2:
nn+1
u=1.0;
j=1;
while（j<
k）
u=u*（x0（i）-x（j））;
j=j+1;
t=t+yy（k）*u;
s（i）=t;
Newton（x,y,x0,3）
1.1619
Newton（x,y,x0,4）
1.1618
（7）实验体会：
本实验通过对牛顿插值公式的程序化，是我们知道了可以用插值公式来解决函数问题，虽然我们只练习了牛顿公式，但是我们对插值公式有了很好的理解，对于在实验过程中遇到的问题都在老师的细心辅导下得到了很好的解决，本次试验受益匪浅。
实验五最小二乘法
5.1实验目的
1掌握最小二乘法的基本思路和拟合步骤；
5.2实验要求
5.3实验内容
（1）题目：
已知一组实验数据如下：
X
2
4
Y
1.95
3.05
3.55
3.85
求问题的最小二乘法。
已知数据对，求多项式
使得
为最小，这就是一个最小二乘问题。
用线性函数为例，拟合给定数据。
算法描述：
步骤1：
输入值，及。
步骤2：
建立法方程组。
步骤3：
解法方程组。
步骤4：
输出。
x=1:
1:
4;
y=[1.95,3.05,3.55,3.85];
p=Polyfit（x,y,3）;
x1=1:
0.1:
y1=polyval（p,x1）;
plot（x,y,'
*r'
x1,y1,'
b'
0.0667-0.70002.7333-0.1500
（6）图形（如果可视化）
通过本次试验，我们发现MATLAB工具箱的功能是很强大的，MATLAB工具箱提供了最小二乘拟合的专门的函数，我们可以通过调用相应的函数就可以达到我们需要的拟合结果，但要了解最小二乘拟合的思想，所以老师建议我们自己编写相应的代码，通过本次试验，也极大的激发了我们的编程兴趣，是我们认识到了我们的不足之处，我们编写的程序都是在正常的运行情况下基本不出错，但是但出现异常时整个程序就会崩溃，也就是说我们编写的程序不够健壮，这需要我们在以后的学习中不断进取。
实验六复化求积公式与龙贝格算法
6.1实验目的
1掌握复化求积公式与龙贝格算法的基本思路和迭代步骤；
6.2实验要求
编写程序，要求给出相应习题的实验结果。
6.3实验内容
用龙贝格算法计算：
龙贝格算法利用外推法，提高了计算精度，加快了收敛速度。
对每一个从2做到，一直做到小于给定的精度是停止计算。
其中（复化梯度求积公式），
输入区间端点，精度控制值，循环次数，定义函数，取，
forto
数据积分近似值。
functions=romberg1（fun,a,b,ep）
%龙贝格算法--求积公式
%此算法只外推到Romberg序列
t1=10000;
t2=-10000;
n=0;
m=1;
h=b-a;
t（1,1）=0.5*（b-a）*（feval（fun,a）+feval（fun,b））;
whileabs（t2-t1）>
=ep
area=0.0;
n=n+1;
h=h/2;
area=area+feval（fun,h*（2*i-1）+a）;
t（n+1,1）=0.5*t（n,1）+area*h;
m=2*m;
ifn>
4
n-j
t（i,j+1）=（4^（j）*t（i+1,j）-t（i,j））/（4^（j）-1）;
t1=t（n-4,4）;
t2=t（n-3,4）;
end
用Romberg序列求得积分近似值为'
s=t2;
functiony=fun（x）
ifx==0
y=1;
y=sin（x）/x
romberg1（'
fun'
e-6）
0.26
以前不知道用程序怎么求积分，但通过本次实验我知道了利用程序是可以求积分的，它主要用的的有限次的迭代，将连续的无限多点的问题有线化，然后就可以利用计算机来解决积分问题，本程序由于比较复杂，所以老师给出了，虽然程序是由老师给出的，但是我们在调用该程序时也基本明白了程序实现的，在实验过程中我们仔细阅读了程序，所以基本上可以再现，但我个人觉得这还远远不够，我们需要自己独立思考学习编写，不断的提高自己的编程能力。
实验七常微分方程数值解法
7.1实验目的
掌握常微分方程数值解的常用算法；
求解常微分方程初值、边值问题的解.
7.2实验题目：
下述方面的相应习题
（1）用改进的欧拉公式，求解常微分方程初值问题的解
（2）用四阶龙格-库塔公式解初值问题
（3）求解常微分方程边值问题的解
7.3实验要求
（1）选择一种计算机语言设计出改进欧拉法和四阶龙格-库塔法方法求解常微分方程初值问题的程序，观察运行结果.
（2）利用Matlab求解常微分方程初值问题
函数dsolve（）用于求解微分方程.Dy表示：
dy/dt（t为缺省的自变量），Dny表示y对t的n阶导数.
Matlab6.5环境下操作如下：
>
y=dsolve（'
Dy=y*y'
'
y（0）=1'
）%求解题目1
Dy=y/t'
y（2.0）=1'
）%求解题目2
（3）ode23，ode45
（1）利用最小二乘法拟合通过改进欧拉法求出微分方程的一系列数值解的近似函数方程.并利用Matlab的绘图功能画出函数的曲线。
（1）题目1：
y’=x^2+y^2,0<
x<
1,
y（0）=0,
h=0.1
题目2：
y’=3*y/（1+x）,0<
1,
y（0）=1,
h=0.2
（4）对应程序1：
functions=Eulerc（fun,x0,xn,y0,h）
%h为步长，改进的求微分方程数值解的Euler公式。
5
n=（xn-x0）/h;
t1=y0+h*feval（fun,x0,y0）;
x0=x0+h;
t2=y0+h*feval（fun,x0,t1）;
y0=（t1+t2）/2;
y（i）=y0;
s=y;
对应程序2：
functions=RK（fun,x0,xn,y0,h）
%h为步长
T=zeros（1,n+1）
升级会员 