山建大数值分析实验报告.docx

上传人:b****8 文档编号:9347527 上传时间:2023-02-04 格式:DOCX 页数:34 大小:40.17KB
下载 相关 举报
山建大数值分析实验报告.docx_第1页
第1页 / 共34页
山建大数值分析实验报告.docx_第2页
第2页 / 共34页
山建大数值分析实验报告.docx_第3页
第3页 / 共34页
山建大数值分析实验报告.docx_第4页
第4页 / 共34页
山建大数值分析实验报告.docx_第5页
第5页 / 共34页
点击查看更多>>
下载资源
资源描述

山建大数值分析实验报告.docx

《山建大数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《山建大数值分析实验报告.docx(34页珍藏版)》请在冰豆网上搜索。

山建大数值分析实验报告.docx

山建大数值分析实验报告

(此文档为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(j

u=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,0

y(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)图形(如果可视化)

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 工作范文 > 演讲主持

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1