matlab课后习题答案第四章

第4章数值运算

习题4及解答

1根据题给的模拟实际测量数据的一组『和)，Q)试用数值差分diff或

数值梯度gradient指令计算yW,然后把)，Q)和V⑺曲线绘制在同

一张图上，观察数值求导的后果。(模拟数据从prob_data40Lmat

获得)

R目的H

・强调：要非常慎用数值导数计算。

•练习mat数据文件中数据的获取。

•实验数据求导的后果

•杷两条曲线绘制在同一图卜.的一种方法。

R解答H

(1)从数据文件获得数据的指令

假如prob_data401.mat文件在当前目录或搜索路径上

clear

loadprob_data401.mat

(2)用diff求导的指令

dt=t(2)-t(l);

yc=diff(y)/dt;%注意yc的长度将比y短1

plot(t,y,*b',t(2:end),yc；r,)

gridon

(3)用gradent求导的指令(图形与上相似)

dt=t(2)-«l);

yc=gradient(7)/dt;

plot(t,y,rb,,t,yc,'r')

gridon

Il说明D

•不到万不得已，不要进行数值求导。

•假若一定要计算数值导数，自变量增量出要取得比原有数据相对误差高1、2个量级以

上。

•求导会使数据中原有为噪声放大。

采用数值计算方法，画出y(x)=k也。在［0,10］区间曲线，并

2JO/

计算y(4.5)o

K提示X

•指定区间内的积分函数可用cumtrapz指令给出。

•＞，(4-5)在计算要求不太高的地方可用find指令算得。

R目的X

•指定区间内的积分函数的数值计算法和cumtrapz指令。

•find指令的应用。

R解答U

dt=le-4;

t=O:dc!O;

t=t+(t==O)*eps;

f=sin(0/t;

s=cumtrapz(f)*dt;

plot(t,s,,LineWidth,,3)

ii=find(t==4.5);

s45=s(iO

s45=

1.6541

3求函数/(%)=*%的数值积分s并请采用符号计算

尝试复算。

R提示X

•数值积分均可尝试。

•符号积分的局限性。

R目的X

•符号积分的局限性。

II解答X

dx=pi/2000;

x=O:dx:pi;

s=trapz(exp(sin(x).A3))*dx

s=

5.1370

符号复算的尝试

symsx

f=exp(sin(x)A3);

ss=int(fpc,0,pi)

Warning:Explicitintegralcouldnotbefound.

>InsvJat58

ss=

int(exp($in(x)A3),x=0..pi)

4用quad求取J二/叫sinxg的数值积分，并保证积分的绝对精度

(tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k4-1));

end

[fobjji]=sort(fobj);%将目标值由小到大排列

tmin=tminGi);%使极小值点做与目标值相应的重新排列

fbbj,tmin

(2)最后确定的最小值点

在N=1,2,…,10的不同分割下,经观察，最后确定出

最小值点是-1.28498111480531

相应目标值是-0.18604801006545

(3)采用作图法近似确定最小值点(另一方法)

(A)在指令窗中运行以下指令：

clear

ft=@(t)sin(5*t).^2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t*cos(2*t);

t=-5:0.001:5;

ff=ft(t);

plot(t,ff)

gridon,shg

(B)经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推

向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据

[tmin2,fobj2]=ginput(l)

tmin2=

-1.28500000993975

fobj2=

-0.18604799369136

出现具有相同数值的刻度区域表明已达最小可分辨状态

(4)符号法求最小值的尝试

symst

fts=sin(5*t)A2*exp(0.06*t*t)-1.5*t*cos(2*t)+1.8*abs(t+0.5);

dfdt=diff(fts,t);%求导函数

tmin=solve(dfdt,t)%求导函数的零点

fobj3=subs(fts,t,tmin)到一个具体的极值点

tmin=

-.60100931947716486053884417850955c-2

fobj3=

.89909908144684551670208797723124

R说明X

•最小值是对整个区间而言的，极小值是对邻域而言的。

•在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。这样可以避免

把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比较中确定的。

•作图法求最小值点，很直观。假若绘图时■,自变曷步长取得足够小，那么所求得的最小

值点有相当好的精度，

•符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。

6设今一3华+2)，⑺=l，M0)=l,竿=0,用数值法和符号法求

atatdt

K目的U

•学习如何把高阶微分方程写成一阶微分方程组。

•ode45解算器的导数函数如何采用匿名函数形式构成。

•如何从Qdc45一组数值解点，求指定自变量对应的函数值。

R解答X

(1)改写高阶微分方程为一阶微分方程组

令MQ)=,()，)，，。)=半，于是据高阶微分方程可写出

at

誓…)

dt

(2)运行以下指令求y(t)为数值解

formatlong

ts=[O,l];

y0=[i；0]；

dydt=®(t,y)[y(2);-2*y(l)+3*y(2)+l];%<4>

%匿名函数写成的ode45所需得导数函数

[tt,yy]=ode45(dydt,ts,yO);

y_05=interpl(tt>yy(:,l),0.5；spline'),%用一维插值求y(0.5)

y-05=

0.78958020790127

(3)符号法求解

symst;

ys=dsolveCD2y-3*Dy+2*y=l,,'y(0)=lJ)y(0)=0,；t')

ys_05=subs(ys,t,sym('0.5'))

ys二

1/2-1/2*cxp(2*t)+cxp(t)

ys.05=

.78958035647060552916850705213780

K说明W

•第条指令中的导数函数也可采用M函数文件表达，具体如下。

functionS=prob_DyDt(t,y)

S=>(2);-2*y(l)+3*y⑵+1];

7已知矩阵人二!11项《8),(1)求该矩阵的“值空间基阵"B;(2)

写出“A的任何列可用基向量线性表出”的验证程序(提示：利

用rref检验)。

g的X

•体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。

•利用rrcf检验两个矩阵能否互为表出。

R解答H

(1)A的值空间的三组不同“基”

A=magic(8);%采用8阶魔方阵作为实验矩阵

[R,ci]=rref(A);

Bl=A(:,d)%直接从A中取基向量

B2=orth(A)%求A值空间的正交基

[VQ]=eig(A);

rv=sum(sum(abs(D))>1000*eps);湘E零特征值数就是矩阵的秩

B3=V(:,l:rv)%取A的非零特征值对应的特征向量作基

B1=

6423

95554

174746

402627

323435

412322

491514

85859

B2=

-0.35360.54010.3536

-0.3536-0.3858-03536

-0.3536-0.2315-0.3536

-0.35360.07720.3536

-0.3536-0.07720.3536

-0.35360.2315-0.3536

-0.35360.3858-0.3536

-0.3536-0.54010.3536

B3=

0.35360.62700.3913

0.3536-0.4815-0.2458

0.3536-0.3361-0.1004

0.35360.1906-0.0451

0.35360.0451-0.1906

0.35360.10040.3361

0.35360.24580.4815

0.3536-0.3913-0.6270

(2)验证A的任何列可月B1线性表出

Bl_A=rref([BlA])%若B1_A矩阵的下5行全为0,

%就表明A可以被B1的3根基向量线性表出

B1_A=

10010011001

01001034-3-47

001001-3-445-7

00000000000

00000000000

00000000000

00000000000

00000000000

B2_A=rref([B2rA])

B2_A=

Columns1through7

1.00(X)00-91.9239-91.9239-91.9239-91.9239

01.0000051.8459-51.8459-51.845951.8459

001.00009.8995-7.0711-4.24261.4142

0000000

0000000

0000000

0000000

0000000

Columns8dirough11

-91.9239-91.9239-91.9239-91.9239

51.8459-51.8459-51.845951.8459

-1.41424.24267.071198995

0000

0000

0000

0000

0000

B3_A=rref([B3rA])

B3_A=

Columns1through7

1.00(X)0091.923991.923991.923991.9239

01.0000042.3447-38.1021-33.859429.6168

001.000012.6462-16.8889-21.131525.3741

0000000

0000000

0000000

0000000

0000000

Columns8through11

91.923991.923991.923991.9239

25.3741-21.1315-16.888912.6462

29.6168-33.8594-38.102142.3447

0000

0000

0000

0000

0000

R说明H

•magic(n)产生魔方阵。魔方阵具有很多特异的性质。就其秩而言，当n为奇数时，该矩

阵满秩；当n为4的倍数时，矩阵的秩总是3；当n为偶数但不是4倍数时，则矩阵的

秩等于(n/2+2)。关于魔方阵的有关历史，请见第6.1.3节。

8已知由MATLAB指令创建的矩阵A=gallery⑸，试对该矩阵进行

特征值分解，并通过验算观察发生的现象。

g的X

•展示特征值分解可能存在的数值问题。

•8ndeig是比较严谨的特征值分解指令。

•Jordan分解的作用。

K解答X

(1)特征值分解

A=gallery(5)

[VJ)]=eig(A);

diagCD)'%为紧凑地显示特征值而写

A=

-911-2163-252

70-69141-4211684

-575575-11493451-13801

3891-38917782-2334593365

1024-10242048-614424572

ans=

Columns1through4

-0.0181-0.0054-0.0171i-0.0054+0.0171i0.0144-0.0104i

Column5

0.0144+0.0104i

(2)验算表明相对误差较大

AE=V*D/V

er_AE=norm(A-AE,,fro')/norm(A,，fro，)%相对F范数

AE二

1.0c+004*

Columns1through4

-0.0009+0.0(X)0i0.0011-0.0(X)0i-O.(X)21+O.(XMX)i0.0063-0.00(X)i

0.0070-O.OOOOi-0.0069+O.OOOOi0.0141-0.0000i-0.0421+O.OOOOi

-0.0575+O.OOOOi0.0575-O.OOOOi-0.1149+O.OOOOi0.3451-O.OOOOi

0.3891-0.00(X)i-0.3891+O.OOOOi0.7781-O.OOOOi-2.3343+O.OOOOi

0.1024-O.OOOOi-0.1024+0.000010.2048-O.OOOOi-0.6144+O.OOOOi

Column5

-0.0252+O.OOOOi

0.1684-().()()(X)i

-1.3800+O.OOOOi

9.3359-0.0001i

2.4570-O.OOOOi

er_/\E=

6.931Oe-OO5

(3)一个更严谨的特征值分解指令

[VcrDc,eigc]=condeig(A)%eigp中的高值时，说明相应的特征值不可信。

Vc=

Columns1through4

-O.OOCX)-0.0000+O.OiMJOi-0.(MX)0-O.OOOOi0.0000+O.OOOOi

0.02060.0207+O.OOOOi0.0207-O.(X)O()i0.0207+O.OOOOi

-0.1397-0.1397+O.OOOOi-0.1397-O.OOOOi-0.1397+O.OOOOi

0.95740.95740.95740.9574

0.25190.2519・O.OOOOi0.2519+O.OOOOi0.2519-O.OOOOi

Column5

0.0000-O.OOOOi

0.0207-O.OOOOi

-0.1397-O.OOOOi

0.9574

0.2519+O.OOOOi

Dc=

Columns1through4

-0.0181000

0-0.0054+0.0171i00

00-0.0054-0.017110

0000.0144+0.01()4i

0000

Column5

0

0

0

0

0.0144-0.01()4i

eigc=

1.0e+011*

5.2687

5.2313

5.2313

5.1725

5.1724

(4)对A采用Jordan分解并检验

[yj^J]=jordan(A);%求出准确的广义特征值，使A*VJ=VJ*D成立。

QJ

AJ=VJ*DJ/VJ

er_iAJ=notm(A-AJ^£ro,)/norm(A^fror)

DJ二

01000

00100

00010

00001

00000

1.0c+004*

-O.(X)()9O.(X)11-O.(M)21O.(M)63-0.0252

0.0070-0.00690.0141-0.04210.1684

-0.05750.0575-0.11490.3451-1.3801

0.3891-0.38910.7782-2.33459.3365

0.1024-0.10240.2048-0.61442.4572

er_AJ=

2.0500e-011

K说明H

•指令condcig的第3输出量cigc给出的是所谓的“矩阵特征值条件数”。当特征条件数

与%—相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大于

“几何重数”的特征值。此时，实施Jordan分解更适宜。

•顺便指出：借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同

概念。前者描述特征值的问题，后者描述矩阵逆的问题。

•本例矩阵A的特征值条件数很高，表明分解不可信。验算也表明，相对误差较大c

•当对矩阵A进行Jordan分解时，可以看到，A具有5重根。当对Jordan分解进行验算

时,相对误差很小。

9求矩阵Ax=b的解，A为3阶魔方阵，b是(3x1)的全1列向量。

R提示U

•了解magic指令

•rref用于方程求解。

•矩阵除法和逆阵法解方程。

II目的X

•满秩方阵求解的一般过程。

•rref用于方程求解。

•矩阵除法和逆阵法解方程。

II解答H

A=magic(3);%产生3阶魔方阵

b=ones(3,l);%(3*1)全1列向量

[R,q=rref([A,b])%GaussJordan消去法解方程，同时判断解的唯一性

x=A\b%矩阵除解方程

xx=inv(A)*b%逆阵法解方程

R=

1.0000000.0657

01.000000.0657

001.00000.0657

123

0.0667

0.0667

0.0667

xx=

0.0667

0.0667

0.0667

R说明1

•rrcf指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的3根“坐标向昆”

和(或)C3指示的3根基向量，可见A3满秩，因此方程解唯一。

•在本例情况下，矩阵除、逆阵法、rref法所得解一致。

10求矩阵Ax=b的解，A为4阶魔方阵，b是(4X1)的全1列向量。

K提示X

•用rreF可观察A不满秩，但b在A的值空间中，这类方程用无数解。

•矩阵除法能正确求得这类方程的特解。

•逆阵法不能求得这类方程的特解。

•注意特解和齐次解

R目的X

•A不满秩，但b在A的值空间中，这类方程的求解过程。

•rref用于方程求解。

•矩阵除法能正确求得这类方程的特解。

•逆阵法不能求得这类方程的特解。

•解的验证方法。

•齐次解的获取。

•全解的获得。

R解答X

(1)借助增广矩阵用指令rref求解

A=ma^c(4);%产生3阶魔方阵

b=ones(4,l)；%全1列向量

[R,q=rref([A,b])%求解，并判断解的唯一性

R=

1.0000001.00900.0588

01.000003.00300.1176

001.0000-3.0000-0.0588

00000

c=

123

关于以上结果的说明:

•R阶梯阵提供的信息

■前4列是原A阵经消元变换后的阶梯阵；而第5列是原b向量经相同变换后的结

果。

■R的前三列为“基"，说明原A阵秩为3;而笫4列的前三个元素，表示原A阵的

第4列由其前三列线性组合而成时的加权系数，即方程的一个解。

■R的第5列表明：b可由原A阵的前三列线性表出；b给出了方程的一个解；由于

原A阵“缺秩”，所以方程的确解不唯一。

•C数组提供的信息

■该数组中的三个元素表示变换取原A阵的第1,2,3列为基。

■该数组的元素总数就是“原A阵的秩”

(2)矩阵除求得的解

x=A\b

Warning:Matrixisclosetosingularorbadlyscaled.

Resultsmav■beinaccurate.RCOND=1.306145e-017.

X=

0.0588

0.1176

-0.0588

0

运行结果指示：矩阵除法给出的解与rref解相同°(实际上，MATLAB在设计“除法”

程序时，针对“b在A值空间中I”的情况，就是用rref求解的。)

(3)逆阵法的解

xx=inv(A)*b

Warning:Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=l.306!45c-017.

xx=

0.0469

0.1875

-0.0625

-0.0156

(3)验证前面所得的解

b_rref=A(:,C)*R(C,5)%验算rref的解

b_d=A*x%验算矩阵除的解

b_inv=A*xx%验算逆阵法的解

b_rref=

1

1

1

1

b_d=

1

1

1

1

b_inv=

0.7344

1.5469

1.1719

1.8594

显然，在本例中，逆阵法的解是错误的。原因是：A不满秩，A的逆阵在理论上不存在。

这里所给出的逆阵是不可信的。

(4)求齐次解

xg=nuU(A)%Ax=O的齐次解

xg=

0.2236

0.6708

-0.6708

-0.2236

(5)方程的全解

齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。

mgdefault%为本书结果可被读者核定而设，并非必要。

f=randn(l,6)%6个随机系数

xx=repmat(x,l,6)+xg*f%产生6个不同的特解

A*xx%所得结果的每列都应该是全1,即等于b.

f=

0.53771.8339-2.25880.86220.3188-1.3077

XX=

0.17900.4689-0.44630.25160.1301-0.2336

0.47831.3479-1.39760.69600.3315-0.7596

-0.4195-1.28901.4565-0.6372-0.27270.8184

-0.1202-0.41010.5051-0.1928-0.07130.2924

ans=

1.00001.00001.00001.00001.00001.0000

1.00001.00001.00001.00001.00001.0000

1.00001.00001.00001.00001.00001.0000

l.OO(X)l.OO(X)l.(X)()()1.0000l.(X)()()1.(X)()()

K解答X

•(在用除法和逆阵法求解时出现)警告信息中RCQND=L306145e-017是矩阵A的估

计条件倒数。该数愈接近0,A就愈“病态”；该数接近1时，A就愈“良态”。该条

件数由rcond(A)给出c注意：rcond条件倒数与cond条件数的算法不同。

1

9

11求矩阵Ax=b的解，A为4阶魔方阵，b=:o

4

II提示X

•由rrcf可以看出A不满秩，b不在A的值空间中，方程没有准确解。

•但可求最小二乘近似解。

咐的X

•A不满秩，b不在A的值空间中，方程没有准确解。

R解答H

(1)借助增广矩阵用指令rrcf求解

A=magjc(4);%产生3阶魔方阵

b=(l：4),;

[R,Q=rref([A,b])%求解，并判断解的唯一性

R=

10010

01030

001-30

00001

c=

1235

(2)用伪逆求最小二乘近似解(超出范围，仅供参考。)

x=pinv(A)*b%非准确解

b_pinv=A*x%验算

0.0235

0.1235

0.1235

0.0235

b_pinv=

1.3000

2.9000

2.10(H)

3.7000

R解答X

•C表明，A的秩为3,A不满秩；R第5列第4元素非零，说明b不在A的值空间中,

因此方程没有准确解，但可以求最小二乘近似解。

12求一0.5+/-1(k^lsinfsinzll=0的实数解。

R提示X

•在适当范围内，作图观察一元复杂函数的形态：观察解的存在性；解的唯一性。进而,

借助图形法求近似解，

•匿名函数的使用方法，

•fzero指令的用法。

K目的X

•作图法求一元复杂函数解上的作用：观察解的存在性；解的唯一性；得近似解。

•匿名函数的使用方法，

•fzero指令的用法。

R解答H

(1)作图观察函数并求近似解

t=-l:0.001:5;

y=@(t)-0.5+t-10*exp(-0.2*t).*abs(sin(sin(t)));

plot(t,y(t))%利用匿名函数求y函数值

gridon,shg

[ttl,yyl]=ginput(l)%从图形获得近似解

ttl=

2.7370

yyi=

0.0097

(2)进一步利用fzero求精确解

[t,yt]=fzero(y,ttl)

2.7341

yt=

2.2204e-015

K说明X

•假如在从图上获取数据前，先把零点附近图形放大：可以得到精度更高的近似解

吗的解。

13求解二元函数方程组f

cosU+y)=0

R目的X

•solve指令的用法。

R解答X

(1)符号法只能得到两组解

S=solveCsin(x-y)'；cos(x+y),；x,；y,);

X=S.x,Y=S.y

X=

[l/4*pi]

Y=

[l/4*pi]

[-l/4*pi]

(2)数值法可以看到许多解，并从中找到规律

x=-6^).l:6;y=x；pC,Y]=meshgnd(x,y);

Zl=sin(X-Y);

Z2=cos(X+Y);

contour(X,Y^l,3)%在Z1的取值范围内取3个等分点作为等位线取值。

%由于Z1关于zl=O对称，所以中间线，即Z1=O的等位线。

axissquare

holdon

contour(X,Y^2,3)%注意：所有绿线交点都是解

holdoff

gridon;shg

14假定某窑工艺瓷器的烧制成品合格率为0.157,现该窑烧制100

件瓷器，请画出合格产品数的概率分布曲线。

R目的JI

•指令binopdf的用法。

•绘图指令stem的用法。

II解答H

N=100;p=0.157;%给定二项分布的特征参数

k=O:N;%定义事件A发生的次数数组

pdf=binopdf(k»N,p);%算出各发生次数的概率

stem(k,pd^

gridon

shg

15试产生均值为4,标准差为2的(10000X1)的正态分布随机数组a,

分别用hist和histfit绘制该数组的频数直方图，观察两张图形的差

异。除histfit上的拟合红线外，你能使这两个指令汇出相同的频

数直方图吗？

K目的X

•指令normrnd的应用，及随机发生器的初始状态控制。

•指令hist与histfit的异同。

R解答X

(1)绘制频数直方图

mgdefault%为本例结果可被读者重现而设，并非必要。

mu=4;sigma=2;%定义均值和标准差

%a=4+2*randn(10000,：);

a=nonnmd(mu,sigma,10000,1)；%产生正态分布随机数组a<3>

subplot(2,l,l)^ist(a)腌制a的频数统计直方图

subplot(2,l,2),histfit(a)

(2)对两个绘图指令中的直方条的数目设置相同的值

subplot(2,l,l)Jiist(a,50)%绘制a的频数统计直方图

subplot(2,l»2)»histfit(a,50)

II说明X

•第V3＞条指令可用一下指令等价代替

a=4+2*randn(10000,l);

16从数据文件prob_data416.mat得到随机数组R,下面有一段求取

随机数组全部数据最大值、均值和标准差的程序。

Mx=max(max(R))JM[e=mean(mean(R)),St=std(std(R)),

试问该程序所得的结果都正确吗？假如不正确，请写出正确的程

序。

R目的U

•mat数据文件中数据的获取。

•max,mean,std的正确运用。

R解答H

(1)获取随机数组R

先要保证数据文件pmb_data416.mat在当前目录上或搜索路径上，然后运行以下指令。

loadprob_data416

(2)运行题给程序

Mx=max(max(R))

Me=mean(mean(R))

S=std(std(R))

Mx=

0.9997

Me=

0.5052

S=

0.0142

(3)正确程序

Mxl=max(R(:))

Mel=mean(R(:))

Sl=std(R(:))%<7>

Mxl=

0.9997

Mel=

0.5052

Si二

0.2919

R说明》

・两种方法得到的最大值和均的计算结果相同。求取R全部数据标准差的正确方法见第

条指令。

•题给方法求出的S不是R全部数据的标准差，而是对于“各列数据样本标准差”的标

准差。该S很小，说明各列求得的标准差相互之间离散度很小。

17已知有理分式蚣)=2,其中N(x)=(3/+x)(x3+0.5),

D(x)=(x2+2x-2)(5/+2/+1)。(1)求该分式的商多项式2(A)和

余多项式“X)o(2)用程序验算D(x)Q(x)+r(.r)=Na)是否成立。

R目的D

•指令conv,deconv的应用。

•计算结果的演算。

II解答H

(1)求商及余

clear

N=conv([3010],[1000.5])

D=conv([l2-2],[5201])

[Q^]=deconv(NJD)

N=

3.000001.00001.500000.50000

D=

512-6-32-2

Q=

0.6000-1.4400

r=

".00()00.000021.880()-5.340()-5.52(X)4.5800-2.8800

(2)验算

NNl=conv(QJD);

m=length(r);

NNl(end-m+l：end)=NNl(end-m+l:end)+r%关键指令

err=norm(N-NNl)/norm(N)%多项式系数向量的相对范数误差

NN1=

3.000001.00001.500000.50000

0

18现有一组实验数据x,y,试求这组数据的拟合多项式。

R目的D

•prob_data418.mat数据文件中数据的获取。

•多项式拟合指令polyfit和多项式求值指令polp-al的应用。

R解答X

loadprob_data418%要注意搜索路径

P=polyRt(玛y,5)%采用5阶多项式拟合

xx=linspace(x⑴Rend),5*length(x));%或为应用需要，或为精细绘图，增加数据点

yy=polyval(P»xx);

plot(x,y,k.',xx,yy,'-b')

P=

-0.00390.0338-0.0227-0.44560.9590-0.0364

K说明D

・拟合多项式的阶数要适当，不宜