数学实验报告2.docx

数学实验报告2.docx

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

数学实验报告2

数学实验报告

院系计算机学院

班级2004级二学位

指导教师张兴永

姓名田琳

实验目录

1．计算

程序内容：

源程序：

\s.m；

程序清单：

functions（x）

s=1;

for（i=1:

8）;

s=s*cos（2^i*x）;

end

s

实验结果：

S

（1）=3.4160e-004;s

（2）=8.1827e-004;s（pi）=1;s（pi/2）=-1;s（pi/3）=0.0039

2.求

的fourier、Laplace、ztrans变换

程序内容：

fourier、Laplace、ztrans变换

源程序：

\bianhuan.m；

程序清单：

symst;

f=5*sin（2*t）-3*cos（2*t）

fourier_f=fourier（f）

laplace_f=laplace（f）

ztrans_f=ztrans（f）

实验结果：

f=5*sin（2*t）-3*cos（2*t）

fourier_f=pi*（5*i*Dirac（w+2）-3*Dirac（w+2）-5*i*Dirac（w-2）-3*Dirac（w-2））

laplace_f=10/（s^2+4）-3*s/（s^2+4）

ztrans_f=

10*z*cos

（1）*sin

（1）/（-4*z*cos

（1）^2+z^2+2*z+1）-3*（z+1-2*cos

（1）^2）*z/（-4*z*cos

（1）^2+z^2+2*z+1）

3．矩阵高斯消去法

对矩阵

进行高斯消去法变换

求三次初等变换的总的等价乘子,并用MATLAB求原来A的行列式、秩和迹

源程序：

\gaosi.m

程序清单：

A=[107;415;2-19];A0=A;%输入A，并保留一个备份

A（2,:

）=-4*A（1,:

）+A（2,:

）

A1=A,

B1=A1/A0%消去A（2,1）,求B1

A（3,:

）=-2*A（1,:

）+A（3,:

）

A2=A,

B2=A2/A1%消去A（3,1）

A（3,:

）=-A（3,2）/A（2,2）*A（2,:

）+A（3,:

）

A3=A,

B3=A3/A2%消去A（3,2）

B0=A3/A0%求三次初等变换的总的等价乘子

det_A=det（A0）,

rank_a=rank（A0）,

tr_A=trace（A0）,%求原来A的行列式、秩和迹

实验结果：

，

det_A=-28rank_a=3tr_A=11

4．画出心型线、星型线、双纽线和四叶玫瑰线的图形

心形线绘制：

Clear

t=0:

0.001:

2*pi;

subplot（2,2,1）;

polar（a,1+cos（t））

subplot（2,2,2）;

plot（cos（t）.^3,sin（t）.^3）

subplot（2,2,3）;

polar（t,abs（sin（t）.*cos（t）））

subplot（2,2,4）;

polar（t,（cos（2*t））.^0.5）

图像：

图表1四曲线

5．求fibonacci数

图表2fib曲线

实验结果：

100以内的fibonacci数：

f=

1123581321345589

实验结果：

1000以内的fibonacci数：

f=

1123581321345589

144233377610987

6．数据三次拟合曲线

数据：

1234567899123165198243277353345303288275的三次拟合

程序：

\nihe3.m

图表3三次拟合曲线

7．递推公式的稳定性

实验内容：

[1]P11页试验课题1

实验程序:

clear;

clc;

symsx;

resualt=zeros（4）;

N_R=1;

a=input（'输入a的值:

'）;

%------------------------------------方案1

I1=log（a+1）-log（a）;

forn=1:

10

I1=-1*a*I1+1/n;

f=x^n/（a+x）;

I0=int（f,'x',0,1）;

I0=vpa（I0,500）;

I0=vpa（I0,8）;

I=vpa（I1,6）;

wucha1=abs（（I0-I1）/I0）;

wucha1=vpa（wucha1,8）;

db=[nI0Iwucha1];

ifn==1

resualt=db;

else

resualt=[resualt;db];

end

end

resualt

%------------------------------------方案2

resualt=zeros（4）;

N=13;

ifa>=N/（N+1）

I2=（2*a+1）/（2*a*（a+1）*（N+1））;

else

I2=0.5*（1/（（a+1）*（N+1））+1/N）;

end

forn=N:

-1:

1

f=x^n/（a+x）;

I0=int（f,'x',0,1）;

I0=vpa（I0,500）;

I0=vpa（I0,8）;

I2=（-1*I2+1/n）/a;

I=vpa（I2,6）;

wucha2=abs（（I0-I2）/I0）;

wucha2=vpa（wucha2,8）;

db=[nI0Iwucha2];

ifn==N

resualt=db;

else

resualt=[resualt;db];

end

end

resualt

实验结果:

输入a的值:

0.05

次数精确值迭代值相对误差

[1,.84777388,.847774,.22248517e-8]

[2,.45761131,.457611,.85349539e-8]

[3,.31045277,.310453,.63500227e-8]

[4,.23447736,.234477,.68175840e-8]

[5,.18827613,.188276,.10198168e-7]

[6,.15725286,.157253,.44935959e-9]

[7,.13499450,.134994,.10844168e-8]

[8,.11825028,.118250,.42221299e-7]

[9,.10519860,.105199,.25088308e-7]

[10,.94740070e-1,.947401e-1,.13928927e-8]

次数精确值迭代值相对误差

[13,.72971840e-1,.889587e-1,.21908205]

[12,.79024729e-1,-.112507,2.4236878]

[11,.86172087e-1,4.06831,46.211490]

[10,.94740070e-1,-79.3663,838.72635]

[9,.10519860,1589.55,15108.966]

[8,.11825028,-31788.4,268824.43]

[7,.13499450,635772.,4709611.4]

[6,.15725286,-.127154e8,80859783.]

[5,.18827613,.254309e9,.13507216e10]

[4,.23447736,-.508617e10,.21691531e11]

[3,.31045277,.101723e12,.32766162e12]

[2,.45761131,-.203447e13,.44458454e13]

[1,.84777388,.406894e14,.47995561e14]

输入a的值:

15

次数精确值迭代值相对误差

[1,.31922183e-1,.319222e-1,.19912866e-8]

[2,.21167256e-1,.211673e-1,.21971093e-8]

[3,.15824494e-1,.158245e-1,.19548773e-8]

[4,.12632590e-1,.126326e-1,.36732305e-7]

[5,.10511157e-1,.105112e-1,.37710110e-8]

[6,.89993133e-2,.899931e-2,.11542747e-6]

[7,.78674434e-2,.786746e-2,.19750588e-5]

[8,.69883483e-2,.698812e-2,.33252485e-4]

[9,.62858867e-2,.628937e-2,.55451370e-3]

[10,.57116994e-2,.565942e-2,.91538526e-2]

次数精确值迭代值相对误差

[13,.44830339e-2,.482067e-2,.75313181e-1]

[12,.48293362e-2,.523418e-2,.83829670e-1]

[11,.52335998e-2,.571166e-2,.91344598e-1]

[10,.57116994e-2,.628589e-2,.10052873]

[9,.62858867e-2,.698835e-2,.11175216]

[8,.69883483e-2,.786744e-2,.12579441]

[7,.78674434e-2,.899931e-2,.14386756]

[6,.89993133e-2,.105112e-1,.16799544]

[5,.10511157e-1,.126326e-1,.20182674]

[4,.12632590e-1,.158245e-1,.25267218]

[3,.15824494e-1,.211673e-1,.33762608]

[2,.21167256e-1,.319222e-1,.50809264]

[1,.31922183e-1,.645385e-1,1.0217452]

明显,对方案2计算结果都不可靠.

8．迭代法的收敛性与收敛速度的比较

实验内容：

[1]P37页试验课题二

实验程序

clear;

clc;

symsx;

f=x^3-sin（x）-12*x+1;

resualt=solve（'x^3-sin（x）-12*x+1','x'）;

disp（'Matlab求根结果:

'）

r=vpa（resualt,8）

pause;

df=diff（f,'x'）;

x0=input（'输入迭代初值:

'）;

N=input（'输入最多迭代次数:

'）;

e=input（'输入迭代精度:

'）;

fork=1:

N

f0=subs（f,x0）;

df0=subs（df,x0）;

ifdf0==0

disp（'导数为0,停止计算'）

break;

else

xx=x0-f0/df0;

ifabs（xx）<=1

E=abs（xx-x0）;

else

E=abs（（x0-xx）/xx）;

end

ifE

disp（'牛顿法__计算结果:

'）

x=vpa（xx,8）

disp（'计算误差为:

'）

error=vpa（abs（（xx-r）/r）,5）

f_x=vpa（subs（f,xx）,8）

disp（'最终迭代次数:

'）

k

break;

else

x0=xx;

end

end

end

ifk==N

disp（'经过设置的迭代次数没有收敛,计算失败'）

end

pause;

disp（'普通迭代法:

'）

disp（'迭代式1:

'）

x0=rand（-4,-3）;

f1=（12*x+sin（x）-1）^（1/3）;

fork=1:

N

xx=subs（f1,x0）;

ifabs（xx）<=1

E=abs（xx-x0）;

else

E=abs（（x0-xx）/xx）;

end

ifE

disp（'[-4,-3]计算结果:

'）

x=vpa（xx,8）

disp（'计算误差为:

'）

error=vpa（abs（（xx-r）/r）,5）

f_x=vpa（subs（f,xx）,8）

disp（'最终迭代次数:

'）

k

break;

else

x0=xx;

end

end

ifk==N

disp（'经过设置的迭代次数没有收敛,计算失败'）

end

pause;

x0=rand（3,4）;

f1=（12*x+sin（x）-1）^（1/3）;

fork=1:

N

xx=subs（f1,x0）;

ifabs（xx）<=1

E=abs（xx-x0）;

else

E=abs（（x0-xx）/xx）;

end

ifE

disp（'[3,4]牛顿法__计算结果:

'）

x=vpa（xx,8）

disp（'计算误差为:

'）

error=vpa（abs（（xx-r）/r）,5）

f_x=vpa（subs（f,xx）,8）

disp（'最终迭代次数:

'）

k

break;

else

x0=xx;

end

end

ifk==N

disp（'经过设置的迭代次数没有收敛,计算失败'）

end

pause;

disp（'普通迭代法:

'）

disp（'迭代式2:

'）

x0=rand（0,0.2）;

f1=（12*x+sin（x）-1）^（1/3）;

fork=1:

N

xx=subs（f1,x0）;

ifabs（xx）<=1

E=abs（xx-x0）;

else

E=abs（（x0-xx）/xx）;

end

ifE

disp（'[0,0.2]计算结果:

'）

x=vpa（xx,8）

disp（'计算误差为:

'）

error=vpa（abs（（xx-r）/r）,5）

f_x=vpa（subs（f,xx）,8）

disp（'最终迭代次数:

'）

k

break;

else

x0=xx;

end

end

ifk==N

disp（'经过设置的迭代次数没有收敛,计算失败'）

end

实验结果

函数

的图像为：

图表4

方案

初值

迭代次数

迭代结果

误差

牛顿法

-3

5

-3.4911788

0

普通迭代

（1）

-3

发散

普通迭代

（2）

-3

发散

牛顿法

3

5

3.4101251

0

普通迭代

（1）

3

发散

普通迭代

（2）

3

发散

牛顿法

-1

4

.76963989e-1

.74718010e-13

牛顿法

0

4

.76963989e-1

0

普通迭代

（1）

0

2

.76964248e-1

-.33644973e-5

[34]

0

2

.77153212e-1

-.24559693e-2

普通迭代

（2）

0

2

.14285360

-.85369616

图表5

9．雅可比迭代法与高斯塞德尔法的收敛性与收敛速度

实验内容：

[1]P71页试验课题（四），雅可比迭代法与高斯——塞德尔迭代法的收敛性与收敛速度。

实验程序：

clear;

clc;

A=input（'输入系数矩阵A:

'）;

b1=input（'输入矩阵b1:

'）;

b2=input（'输入矩阵b2:

'）;

disp（'Matlab计算结果:

'）

x1=inv（A）*b1;

X1=vpa（x1,6）

x2=inv（A）*b2;

X2=vpa（x2,6）

L=tril（A）;

U0=triu（A）;

U=L-A;

L=U0-A;

D=A+L+U;

pause;

N=input（'输入最多迭代次数:

'）;

disp（'雅可比迭代计算结果:

'）

MJ=inv（D）*（L+U）;

disp（'计算b1情况'）;

X0=zeros（sqrt（numel（A））,1）;

%与输入矩阵配置相同的初始值

fork=1:

N

X=MJ*X0+D^-1*b1;

if（norm（X-X0）<1e-6）

b1

disp（'计算结果:

'）

X1=vpa（X,6）

disp（'迭代次数'）

k

AX=vpa（A*X,6）

break;

else

X0=X;

end

end

ifk==N

disp（'迭代没有收敛'）

end

pause;

disp（'计算b2情况'）;

X0=zeros（sqrt（numel（A））,1）;

%与输入矩阵配置相同的初始值

fork=1:

N

X=MJ*X0+D^-1*b2;

if（norm（X-X0）<1e-6）

b2

disp（'计算结果:

'）

X2=vpa（X,6）

disp（'迭代次数'）

k

AX=vpa（A*X,6）

break;

else

X0=X;

end

end

ifk==N

disp（'迭代没有收敛'）

end

pause

disp（'高斯--塞德尔迭代计算结果:

'）

MG=（D-L）^-1*U;

disp（'计算b1情况'）;

X0=zeros（sqrt（numel（A））,1）;

%与输入矩阵配置相同的初始值

fork=1:

N

X=MG*X0+（D-L）^-1*b1;

if（norm（X-X0）<1e-6）

b1

disp（'计算结果:

'）

X1=vpa（X,6）

disp（'迭代次数'）

k

AX=vpa（A*X,6）

break;

else

X0=X;

end

end

ifk==N

disp（'迭代没有收敛'）

end

pause;

disp（'计算b2情况'）;

X0=zeros（sqrt（numel（A））,1）;

%与输入矩阵配置相同的初始值

fork=1:

N

X=MG*X0+（D-L）^-1*b2;

if（norm（X-X0）<1e-6）

b2

disp（'计算结果:

'）

X2=vpa（X,6）

disp（'迭代次数'）

k

AX=vpa（A*X,6）

break;

else

X0=X;

end

end

ifk==N

disp（'迭代没有收敛'）

end

实验结果

1．

，

，

结果：

，

雅可比迭代法计算b1迭代24次，计算b2迭代30次。

高斯塞德尔法计算b1迭代15次，计算b2迭代20次。

2．

，

，

结果：

雅可比迭代法计算计算不收敛。

高斯塞德尔法计算b1迭代45次，计算b2迭代52次。

3．

，

结果为

但使用两种方法都不收敛。

10．龙格现象的发生、防止和插值效果的比较

实验内容：

[1]P102页实验课题

（一）龙格现象的发生、防止和插值效果的比较

实验程序：

clear;

clc;

f='x/（1+x^4）';

N=input（'输入插值次数'）;

xi=zeros（1,N+1）;

h=10/N;

fork=1:

N+1

xi（k）=-5+（k-1）*h;

end

yi=subs（f,xi）;

xk=zeros（1,41）;

fork=1:

41

xk（k）=-5+0.25*（k-1）;

end;

%拉格朗日插值算法

ykl=zeros（1,41）;

fori=1:

41

ykl（i）=0;

fork=1:

N+1

t=1;

forj=1:

N+1

ifj~=k

t=（xk（i）-xi（j））/（xi（k）-xi（j））*t;

ykl（i）=ykl（i）+t*yi（k）;

end

end

end

end

subplot（2,2,1）;

plot（xk,ykl,'b'）

title（'拉格朗日插值曲线'）

%分段线性插值

ykx=zeros（1,41）;

fori=1:

41

k=1;

k=k+1;

whilek=xi（k）

k=k+1;

end

ykx（i）=yi（k-1）*（xk（i）-xi（k））/（xi（k-1）-xi（k））+yi（k）*（xk（i）-xi（k-1））/（xi（k）-xi（k-1））;

end

subplot（2,2,2）;

plot（xk,ykx,'g'）

title（'分段线性插值曲线'）

%

yky=zeros（1,41）;

yky=interp1（xi,yi,xk,'spline'）;

subplot（2,2,3）;

plot（xk,yky,'r'）

title（'三次样条插值曲线'）

disp（'计算结果输出:

'）

z=[xk'ykl'ykx'yky']

结算结果

1．

的十次插值图像为：

二十次插值图像：

插值顺序：

，拉格朗日插值、分段线形插值、I型三次样条插值：

十次插值：

-5.0000-0.0799-0.0080-0.0080

-4.75002.7400-0.0099-0.0007

-4.50002.3631-0.0118-0.0011

-4.25001.0023-0.0137-0.0068

-4.0000-0.2309-0.0156-0.0156

-3.7500-0.9262-0.0208-0.0252

-3.5000-1.0791-0.0261-0.0334

-3.2500-0.8821-0.0313-0.0380

-3.0000-0.5846-0.0366-0.0366

-2.7500-0.4067-0.0569-0.0300

-2.5000-0.4943-0.0771-0.0309

-2.2500-0.9042-0.0974-0.0549

-2.0000-1.6096-0.1176-0.1176

-1.7500-2.5185-0.2132-0.2259

-1.5000-3.4987-0