当前位置:文档之家 > matlab实验报告

matlab实验报告

MATLAB 实验报告

1、在区间[-1,1]上分别取n=10、20用两组选中节点对龙格函数22511)(x

x f +=

作多项式插值及三次样条插值,对每个n 值,分别画出插值函数及f(x)的图形.

解:

n=10

在matlab 命令窗口中键入:

>>x=-1:0.2:1;

y=1./(1+25*x.^2);

y1=interp1(x,y,'pchip');

y2=interp1(x,y,'spline');

plot(x,y1,'o',x,y1,'-',x,y2,'*',x,y2,'-.'),grid

legend('样本点','三次插值','三次样条插值')

回车得出:

matlab实验报告

n=20

在matlab 命令窗口中键入:

>> x=-1:0.1:1;

y=1./(1+25*x.^2);

y1=interp1(x,y,'pchip');

y2=interp1(x,y,'spline');

plot(x,y1,'o',x,y1,'-',x,y2,'*',x,y2,'-.'),grid

legend('样本点','三次插值','三次样条插值')legend('样本点','三次插值','三次样条插值') 回车得出:

matlab实验报告

由结果可见,用两种方法画出的曲线在样本点之间取值并无太大差异,曲线亦基本上一致。

2、对于给定函数22511)(x

x f +=在区间[-1,1]上取)10,,1,0(2.01 =+-=i i x i ,试求3次曲线拟合,试画出拟合曲线并打印出方程。

解:在matlab 命令窗口中键入:

>> x=-1:0.2:1;y=1./(25*x.^2+1);

p=polyfit(x,y,3)

回车得出:

p =

-0.0000 -0.5752 0.0000 0.4841

即拟合的多项式为:4814.05752.02

+-=x y

键入:

x1=-1:0.1:1;y1=polyval(p,x1);

plot(x,y,'*',x1,y1),grid

legend('样本点','拟合曲线p3(x)')

回车得出:

matlab实验报告

由结果可看到,拟合曲线并未通过所有的样本点,它冲出了一些随机误差,更能真实地反映出两组量间关系变化的趋势。

3、用LU 分解及列主元高斯消去法解线性方程组

?????

???????=????????????????????????15900001.582 0 1 21- 5 1- 52 6 2.099999 3-1 0 7- 104321x x x x . 输入Ax=b 中系数A=LU 分解的矩阵L 及U,解向量x 及detA;列主元法的行交换次序,解向量x 及detA;比较两种方法所得结果。

解:LU 分解法

在matlab 命令窗口中键入:

>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];

b=[8 5.900001 5 1]';

[L,U]=lu(A)

x=U\(L\b)

回车得出:

L =

1.0000 0 0 0

-0.3000 -0.0000 1.0000 0

0.5000 1.0000 0 0

0.2000 0.9600 -0.8000 1.0000

U =

10.0000 -7.0000 0 1.0000

0 2.5000 5.0000 -1.5000

0 0 6.0000 2.3000

0 0 0 5.0800

x =

0.0000

-1.0000

1.0000

1.0000

即解向量x=[0 -1 1 1]’

键入:

>>det(A)

回车得出:

ans =

-762.0001

即行列式det(A)=-762.0001.

列主元高斯消去法

在matlab命令窗口中键入:

>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];

b=[8 5.900001 5 1]';

x=gauss(A,b)

回车得出:

x =

0.0000 -1.0000 1.0000 1.0000

即解向量x=[0 -1 1 1]’

键入:

>> det(A)

回车得出:

ans =

-762.0001

即行列式det(A)=-762.0001.

由结果可见,用LU分解法与列主元消去法所得到的结果一致。

4、 线性方程组Ax=b 的A 及b 为

?????

???????=10 9 5 7 9 10 6 8 5 6 5 7 7 8 7 10A ,????????????=31332332b , 则解().1,1,1,1T

x =用MATLAB 内部函数求detA 及A 的所有特征值和cond(A)2. 解:在matlab 命令窗口中键入:

>>A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];

det(A)

[V,L]=eig(A)

c=cond(A)

回车得出:

ans =

1

V =

0.5016 -0.3017 0.6149 0.5286

-0.8304 0.0933 0.3963 0.3803

0.2086 0.7603 -0.2716 0.5520

-0.1237 -0.5676 -0.6254 0.5209

L =

0.0102 0 0 0

0 0.8431 0 0

0 0 3.8581 0

0 0 0 30.2887

c =

2.9841e+003

由结果可知,det(A)=1,A 的四个特征值为0.0102,0.8431,

3.8581,30.2887。cond(A)2=2.9841e+003。

5、给出线性方程组b x H n =,其中系数矩阵n H 为希尔伯特矩阵:

n n ij n R h H ?∈=)(,1

1-+=j i h ij , .,,2,1,,n j i = 假设**,)1,,1,1(x H b R x n n T =∈= ,若取n=6,8,10,分别用雅可比迭代法及SOR 迭代(ω=1,1.25,1.5)求解.比较计算结果.

解:n=6

在matlab 命令窗口中键入:

>>H=hilb(6);

x0=[1 1 1 1 1 1]';

b=H*x0

回车得出:

b =

2.4500

1.5929

1.2179

0.9956

0.8456

0.7365

jacobi 迭代法:

>>H=hilb(6);

b=[2.4500 1.5929 1.2179 0.9956 0.8456 0.7365]';

x0=[1 1 1 1 1 1]';

[x,n]=jacobi(H,b,x0)

回车得出:

Warning:迭代次数太多,可能不收敛!

x =

1.0e+122 *

0.3290

0.7099

0.9462

1.1118

1.2354

1.3317

n =

200

SOR 迭代法:(ω=1)

>> [x,n]=SOR(H,b,x0,1)

回车得出:

Warning:迭代次数太多,可能不收敛!

x =

0.9991

0.9958

1.0392

0.9383

1.0265

1.0001

n =

200

由结果可见,SOR法所给出的解为x=[0.9991,0.9958,1.0392,0.9383,1.0265,1.0001]’,误差比jacobi迭代法要小得多。

n=8

在matlab命令窗口中键入:

>> H=hilb(8);

x0=[1 1 1 1 1 1 1 1]';

b=H*x0

回车得出:

b =

2.7179

1.8290

1.4290

1.1865

1.0199

0.8968

0.8016

0.7254

jacobi迭代法:

>> H=hilb(8);

x0=[1 1 1 1 1 1 1 1]';

b=[2.7179 1.8290 1.4290 1.1865 1.0199 0.8968 0.8016 0.7254]';

[x,n]=jacobi(H,b,x0)

回车得出:

Warning:迭代次数太多,可能不收敛!

x =

1.0e+151 *

-0.8233

-1.8632

-2.5551

-3.0629

-3.4556

-3.7700

-4.0281

-4.2441

n =

200

SOR迭代法:(ω=1.25)

>>[x,n]=SOR(H,b,x0,1.25)

回车得出:

Warning:迭代次数太多,可能不收敛!

x =

0.9982

1.0124

1.0177

0.8680

1.1630

0.8808

1.1393

0.9211

n =

200

由结果可见,SOR法所给出的解为

x=[0.9982,1.0124,1.0177,0.8680,1.1630,0.8808,1.1393,0.8211]’,误差比jacobi迭代法要小得多。

n=10

>> H=hilb(10);

x0=[1 1 1 1 1 1 1 1 1 1]';

b=H*x0

回车得出

b =

2.9290

2.0199

1.6032

1.3468

1.1682

1.0349

0.9307

0.8467

0.7773

0.7188

jacobi迭代法:

在matlab窗口中键入:

>> H=hilb(10);

x0=[1 1 1 1 1 1 1 1 1 1]';

b=[2.9290 2.0199 1.6032 1.3468 1.1682 1.0349 0.9307 0.8467 0.7773 0.7188]';

[x,n]=jacobi(H,b,x0)

回车得出:

x =

1.0e+173 *

-0.2452

-0.5736

-0.8029

-0.9770

-1.1151

-1.2280

-1.3223

-1.4024

-1.4715

-1.5316

n =

200

SOR迭代法:(ω=1.5)

>> H=[x,n]=SOR(H,b,x0,1.5)

回车得出:

Warning:迭代次数太多,可能不收敛!

x =

0.9977

1.0191

0.9711

1.0367

0.8891

1.1433

0.7735

1.2451

1.0809

0.8423

n =

200

由结果可见,SOR法所给出的解为

x=[0.9977,1.0191,0.9711,1.0367,0.8891,1.1433,0.7735,1.2451,1.0809,0.8423]’,误差比jacobi迭代法要小得多。

6、已知矩阵

????????????????????=????????????????=????????????=111 7

1 61 71 31 2161 21 1H , 0 1 0 0 09 8

2 0 08 7 6

3 07 6 5

4 46

5 4 3 2B , 10 9 5 7 9 10

6 8 5 6 5

7 7

8 7 106 A (1) 用MATLAB 函数“eig ”求矩阵全部特征值.

(2) 用基本QR 算法求全部特征值(可用MATLAB 函数“qr ”实现矩阵的QR 分解). 解:(1)

矩阵A

在MATLAB 命令窗口中输入:

>> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];

>> [V,L]=eig(A)

回车得出:

V =

0.5016 -0.3017 0.6149 0.5286

-0.8304 0.0933 0.3963 0.3803

0.2086 0.7603 -0.2716 0.5520

-0.1237 -0.5676 -0.6254 0.5209

L =

0.0102 0 0 0

0 0.8431 0 0

0 0 3.8581 0

0 0 0 30.2887

由计算结果可知,矩阵A 的四个特征值为0.0102,0.8431,3.8581,30.2887.

矩阵B

在MATLAB 命令窗口中输入:

>> B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0];

>> [V,L]=eig(B)

回车得出:

V =

0.4801 0.4841 0.2507 0.4050 -0.2057

0.6623 0.7574 0.5992 -0.7475 0.2129

0.5252 0.3562 -0.7461 0.4928 -0.2385

0.2340 -0.2525 0.1239 0.0673 -0.6295

0.0178 -0.0385 0.0776 -0.1723 0.6776

L =

13.1724 0 0 0 0

0 6.5519 0 0 0

0 0 1.5957 0 0

0 0 0 -0.3908 0

0 0 0 0 -0.9291

由计算结果可知,矩阵B的五个特征值为13.1724,6.5519,1.5957,-0.3908,-0.9291.

矩阵H6

在MATLAB命令窗口中输入:

>> H6=[1 1/2 1/3 1/4 1/5 1/6;1/2 1/3 1/4 1/5 1/6 1/7;1/3 1/4 1/5 1/6 1/7 1/8;1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10;1/6 1/7 1/8 1/9 1/10 1/11];

[V,L]=eig(H6)

V =

-0.0012 -0.0111 0.0622 0.2403 -0.6145 0.7487

0.0356 0.1797 -0.4908 -0.6977 0.2111 0.4407

-0.2407 -0.6042 0.5355 -0.2314 0.3659 0.3207

0.6255 0.4436 0.4170 0.1329 0.3947 0.2543

-0.6898 0.4415 -0.0470 0.3627 0.3882 0.2115

0.2716 -0.4591 -0.5407 0.5028 0.3707 0.1814

L =

0.0000 0 0 0 0 0

0 0.0000 0 0 0 0

0 0 0.0006 0 0 0

0 0 0 0.0163 0 0

0 0 0 0 0.2424 0

0 0 0 0 0 1.6189

由计算结果可知,矩阵H6的六个特征值为0.0000,0.0000,0.0006,0.0163,0.2424,1.6189.

(2)建立M-file,代码如下:

function l = qrtz(A,M)

for (i=1:M)

[q,r]=qr(A);

A=r*q;

l=diag(A);

end

矩阵A

在MATLAB中键入:

>>A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];

l=qrtz(A,20)

回车得出:

l =

30.2887

3.8581

0.8431

0.0102

即矩阵A的四个特征值为30.2887,3.8581,0.8431,0.0102.

矩阵B

在MATLAB中键入:

>> B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0];

l=qrtz(B,20)

回车得出:

l =

13.1723

6.5519

1.5957

-0.9291

-0.3908

即矩阵B的五个特征值为13.1723,6.5519,1.5957,-0.9291,-0.3908.

矩阵H6

在MATLAB中键入:

>> H6=[1 1/2 1/3 1/4 1/5 1/6;1/2 1/3 1/4 1/5 1/6 1/7;1/3 1/4 1/5 1/6 1/7 1/8;1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10;1/6 1/7 1/8 1/9 1/10 1/11];

l=qrtz(H6,20)

回车得出:

l =

1.6189

0.2424

0.0163

0.0006

0.0000

0.0000

即矩阵H6的六个特征值为1.6189,0.2424,0.0163,0.0006,0.0000,0.0000.