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