计算方法实验报告.docx

计算方法实验报告.docx

《计算方法实验报告.docx》由会员分享，可在线阅读，更多相关《计算方法实验报告.docx（18页珍藏版）》请在冰豆网上搜索。

计算方法实验报告

实验一塞德尔迭代法求解线性方程组

【实验目的】

熟悉塞德尔迭代法求解线性方程组的方法

【实验内容】

1.了解matlab语言。

2.用塞德尔迭代法解线性方程组（迭代五次）。

并与精确解

比较。

【实验所使用的仪器设备与软件平台】

计算机、matlab7.0

【实验方法与步骤】

塞德尔迭代算法：

S1：

输入

，允许误差E，最大迭代次数N，迭代初值

；

S2：

计算出

的维数为

；

S3：

置

；

S4：

计算：

S5：

若

则输出

，并结束；

若

，当

时，置

，并转到s4；

当

时，无输出，并结束。

【实验结果】

实验程序：

functionx1=GS（A,b,x0,N,E）

[m,n]=size（A）;%计算A的维数为m

x（m,2）=0;%初始化一个m行2列的矩阵，记录相邻两次迭代结果

x（:

2）=x0;

fori=[1:

1:

N],%迭代最高次数为N

ifmax（abs（x（:

1）-x（:

2）））

break;

else

x（:

1）=x（:

2）;

forj=[1:

1:

n],

ifj==1,

s1=0;

fore1=[2:

1:

m],

s1=s1+A（1,e1）*x（e1,1）;%求和

end

x（j,2）=（b（1,1）-s1）/A（1,1）;

elseifj==m,

s2=0;

fore2=[1:

1:

m-1],

s2=s2+A（m,e2）*x（e2,2）;%求和

end

x（j,2）=（b（m,1）-s2）/A（m,m）;

else

s3=0;

fore3=[1:

1:

j-1],

s3=s3+A（j,e3）*x（e3,2）;%求和

end

fore4=[j+1:

1:

m],

s3=s3+A（j,e4）*x（e4,1）;%求和

end

x（j,2）=（b（j,1）-s3）/A（j,j）;

end

end

end

end

x1=x（:

2）;

【结果分析与讨论】

在matlab命令窗口中输入：

>>A=[5-1-1-1;-110-1-1;-1-15-1;-1-1-110;];

>>b=[-4;12;8;34;];

>>x0=[1;1;1;1];

>>GS（A,b,x0,5,0.000001）

ans=

0.996552514803221

1.998599188612620

2.998420412137192

3.999357211555304

评价：

相对于简单迭代法，赛德尔迭代法收敛速度很快，此程序可以对任意维的方程组求指定的精度的迭代结果，实用性比较强。

实验二最小二乘法拟合数据

【实验目的】

熟悉最小二乘法拟合数据的方法

【实验内容】

1.了解matlab语言

2.给出数据：

-1.00

0.50

-0.50

-0.25

0.00

10.2209

0.3295

0.8826

1.4392

2.0003

0.25

0.50

0.75

1.00

2.5645

3.1334

3.7061

4.2836

用一次、二次多项式利用最小二乘法拟合这些数据，试写出正规方程组，并求出最小平方逼近多项式。

【实验所使用的仪器设备与软件平台】

计算机、matlab7.0

【实验方法与步骤】

最小二乘法拟合算法：

S1：

输入数据列向量

与

，拟合多项式的次数

；

S2：

计算

的列数为

；

S3：

置

（其中

表示矩阵

的第1行第1列，以下的

等类似）；

S4：

计算：

S5：

利用

分解法求解

，得到

次最小平方逼近多项式的系数与多项式；

S6：

在同一坐标系中绘出原数据点与次

多项式逼近曲线。

【实验结果】

实验程序：

functiony=fmintwo（a,m,x）%调用拟合多项式的函数

y=a（1,1）;

fore1=[1:

1:

m],

y=y+（a（e1+1,1））*（x.^e1）;

end

functiona=mintwo（x,y,m）%最小二乘法拟合程序

[n,p]=size（x）;

a（m+1,1）=0;

A（m+1,m+1）=0;

b（m+1,1）=0;

c（2*m+1,1）=0;

c（1,1）=n;

fork=[2:

1:

2*m+1],

s1=0;

fore1=[1:

1:

n],

s1=s1+（x（e1,1））^（k-1）;

end

c（k,1）=s1;

end

fori=[1:

1:

m+1],

forj=[1:

1:

m+1],

A（i,j）=c（i+j-1,1）;%求正规方程组的系数矩阵

end

end

A

forh=[1:

1:

m+1],

s2=0;

fore2=[1:

1:

n],

s2=s2+y（e2,1）*（（x（e2,1））^（h-1））;

end

b（h,1）=s2;%求正规方程组的

向量

end

b

[LDa]=ja2911（A,b）;%调用

分解法求解正规矩阵的程序

z=[min（x）:

0.1:

max（x）];

plot（x,y,'*',z,fmintwo（a,m,z）,'r'）;%绘出原数据表示的点和拟合曲线

附

分解法求系数矩阵为对称矩阵的方程组程序：

function[LDx]=ja2911（A,b）

[n,n]=size（A）;

L=eye（size（A））;

D=zeros（size（A））;

fori=[2:

1:

n],

A（i,1）=A（i,1）/A（1,1）;

forj=[2:

1:

n],

ifi>j;

s1=0;

fore1=[1:

1:

j-1],

s1=s1+A（e1,e1）*A（i,e1）*A（j,e1）;

end

A（i,j）=（A（i,j）-s1）/A（j,j）;

elseifi==j,

s2=0;

fore2=[1:

1:

i-1],

s2=s2+A（e2,e2）*A（i,e2）*A（i,e2）;

end

A（i,j）=A（i,j）-s2;

end

end

end

forh=[1:

1:

n],

D（h,h）=A（h,h）;

end

fork=[2:

1:

n],

fort=[1:

1:

k-1],

L（k,t）=A（k,t）;

end

end

y（n,1）=0;

x（n,1）=0;

y（1,1）=b（1,1）/D（1,1）;

fora=[2:

1:

n],

s3=0;

fore3=[1:

1:

a-1],

s3=s3+L（a,e3）*D（e3,e3）*y（e3,1）;

end

y（a,1）=（b（a,1）-s3）/D（a,a）;

end

x（n,1）=y（n,1）;

forbb=[n-1:

-1:

1],

s4=0;

fore4=[bb+1:

1:

n],

s4=s4+L（e4,bb）*x（e4,1）;

end

x（bb,1）=y（bb,1）-s4;

end

【结果分析与讨论】

在matlab命令窗口中输入：

>>xx=[-1;-0.75;-0.5;-0.25;0;0.25;0.5;0.75;1;];

>>yy=[-0.2209;0.3295;0.8826;1.4392;2.0003;2.5645;3.1334;…

3.7061;4.2836;];

>>mintwo（xx,yy,1）

A=

9.0000000000000000

03.750000000000000

b=

18.118299999999998

8.443674999999999

ans=

2.013144444444444

2.251646666666666

输出的图像：

>>mintwo（xx,yy,2）

A=

9.00000000000000003.750000000000000

03.7500000000000000

3.75000000000000002.765625000000000

b=

18.118299999999998

8.443674999999999

7.586956250000000

ans=

2.000100432900433

2.251646666666666

.0313********

输出的图像：

一次多项式拟合的正规方程组为：

一次拟合多项式为：

二次多项式拟合的正规方程组为：

二次拟合多项式为：

评价：

此程序实现了对数据利用最小二乘法进行

次多项式逼近，并且可以输出原数据表示的点和拟合曲线在同一坐标系中的图像，可以很直观的进行观察拟合的程度。

实验三龙贝格求积分

【实验目的】

熟悉龙贝格求积分的方法

【实验内容】

1.了解matlab语言

2.用龙贝格方法计算积分：

要求误差不超过

。

【实验所使用的仪器设备与软件平台】

计算机、matlab7.0

【实验方法与步骤】

龙贝格求积分算法：

S1：

输入区间左右端点

与

和精度

；

S2：

置

；

S3：

；

S4：

计算：

S5：

置

；

S6：

计算：

S7：

若

，则输出

，并且结束；

若

时，当

则转到s5；

当

则无输出，并且结束。

【实验结果】

实验程序：

functiony=fromberg（x）%可被调用的积分函数

y=2*exp（-（x^2））/（（pi）^（1/2））;

functionI=romberg（a,b,E）%龙贝格求积分程序

T（50,1）=0;

C（50,1）=0;

R（50,1）=0;

h（50,1）=0;

h（1,1）=（b-1）/2;

T（1,1）=h（1,1）*（fromberg（a）+fromberg（b））;

fork1=[2:

1:

5],

h（k1,1）=（b-a）/（2^（k1-1））;

g1=0;

fore1=[1:

1:

2^（k1-2）],

g1=g1+fromberg（a+（2*e1-1）*h（k1,1））;

end

T（k1,1）=（T（k1-1,1）/2）+h（k1,1）*g1;%计算

end

fori1=[1:

1:

4],

S（i1,1）=（4*T（i1+1,1）-T（i1,1））/3;%计算

end

fori2=[1:

1:

3],

C（i2,1）=（16*S（i2+1,1）-S（i2,1））/15;%计算

end

fori3=[1:

1:

2],

R（i3,1）=（64*C（i3+1,1）-C（i3,1））/63;%计算

end

fork=[1:

1:

40],

ifabs（R（k+1,1）-R（k,1））<=E,%误差判断

I=R（k+1,1）;

break;

else%如果不符合精度要求，则自动将区间分半继续计算

h（k+5,1）=（b-a）/（2^（k+4））;

g2=0;

fore2=[1:

1:

2^（k+3）],

g2=g2+fromberg（a+（2*e2-1）*h（k+5,1））;

end

T（k+5,1）=（T（k+4,1）/2）+h（k+5,1）*g2;

S（k+4,1）=（4*T（k+5,1）-T（k+4,1））/3;

C（k+3,1）=（16*S（k+4,1）-S（k+3,1））/15;

R（k+2,1）=（64*C（k+3,1）-C（k+2,1））/63;

end

end

【结果分析与讨论】

在matlab命令窗口中输入：

>>romberg（0,1,10^（-5））

ans=

0.842693582047394

评价：

龙贝格求积分方法精确度比较高，该程序很好的实现了龙贝格方法，但是程序的某些地方还可以简化，使程序精简。

实验四标准四阶R-K方法

【实验目的】

熟悉标准四阶R-K的方法

【实验内容】

1.了解matlab语言

2.用标准四阶R-K求解初值问题

的解在-0.9，-0.8，-0.7的近似值（按四位小数计算）。

【实验所使用的仪器设备与软件平台】

计算机、matlab7.0

【实验方法与步骤】

S1：

输入区间左右端点

与

、

的初值

、初始步长

、最大迭代次数

、允许误差

；

S2：

置

；

S3：

计算：

S4：

置

；

S4：

计算：

S5：

若

，则将

，输出

，并且结束；

若

，当

，将

，并转向s4；

当

，无输出，并且结束。

【实验结果】

实验程序：

functionfk=fRK（x,y）%微分方程的

函数

fk=（x^2）-（y^2）;

functionx=RK（a,b,ya,h）%用R-K方法求微分方程数值解的程序

[m,n]=size（[a:

h:

b]）;

x（n,2）=0;

x（:

1）=[a:

h:

b]';

x（1,2）=ya;

k（4,1）=0;

fori=[2:

1:

n],

k（1,1）=h*fRK（x（i-1,1）,x（i-1,2））;

k（2,1）=h*fRK（x（i-1,1）+h/2,x（i-1,2）+k（1,1）/2）;

k（3,1）=h*fRK（x（i-1,1）+h/2,x（i-1,2）+k（2,1）/2）;

k（4,1）=h*fRK（x（i-1,1）+h,x（i-1,2）+k（3,1））;

x（i,2）=x（i-1,2）+（k（1,1）+2*k（2,1）+2*k（3,1）+k（4,1））/6;

end

function[xh1]=WRK（a,b,ya,h,N,E）%R-K的步长自动选择程序

m（N,1）=0;

x1（N,1）=0;

fori=[1:

1:

N],

[lm（i,1）]=size（[a:

h/（2^（i-1））:

b]）;

end;

y=RK（a,b,ya,h）;

x1（1,1）=y（m（1,1）,2）;

y=RK（a,b,ya,h/2）;

x1（2,1）=y（m（2,1）,2）;

forj=[3:

1:

N],

ifabs（x1（j-1,1）-x1（j-2,1））<=E,

x=RK（a,b,ya,h/（2^（j-2）））;

h1=h/（2^（j-2））;

break;

else

y=RK（a,b,ya,h/（2^（j-1）））;

x（j,1）=y（m（3,1）,2）;

end

end

（该程序输出的是满足精度要求的区间分点的数值解和及相应步长）

【结果分析与讨论】

在matlab命令窗口中输入：

>>[xh]=WRK（-1,-0.5,0,0.1,20,0.00001）

x=

-1.0000000000000000

-0.9500000000000000.047503044462813

-0.9000000000000000.090047822422416

-0.8500000000000000.127736329909335

-0.8000000000000000.160727765024206

-0.7500000000000000.189********6077

-0.7000000000000000.213483856644815

-0.6500000000000000.233766*********

-0.6000000000000000.250369794944545

-0.5500000000000000.263601751827377

-0.5000000000000000.273776792527207

h=

0.050000000000000

所以微分方程的解在-0.9，-0.8，-0.7的近似值分别为0.09004、0.16072、0.21348

评价：

标准四阶R-K方法的精确度比较高，而且算法容易实现。

次程序实现了标准四阶R-K算法，通过多次调用函数的形式使程序得到了精简。