操作3

syms t;
ft=exp(sin(t));
sx=int(ft, t, 0, 4)


delt=2*pi/200;
x=0:delt:2*pi;
y=exp(-sin(x));
sx=delt*cumtrapz(y);
plot(x,y, 'r', 'LineWidth', 6);
hold on;
plot(x, sx, '.b', 'MarkerSize', 15);
plot(x, ones(size(x)), 'k');


*****************************************************************************************
page 3

例3.4.1-1 直接输入向量构造多项式

T=[2,5,0,4,1,4];
fx=poly2sym(T)

****************************************************************************************
page 4
例 4.1.1-2 用多项式的根构造多项式,根为r=[1,2,3,4]
r=[1,2,3,4];
T1=poly(r);
y=poly2sym(T1)
y_class=class(y)

*******************************************************************************************
page 6
例4.1.1-3 多项式的加减乘除运算
f1(x)=2x5+5x4+4x2+x+4, f2(x)=5x3+x2+3x+2

T1=[2,5,0,4,1,4];
T2=[0,0,5,1,3,2];
T3=[5,1,3,2];
T=T1+T2;
T_add=poly2sym(T)

T=T1-T2;
T_sub=poly2sym(T)

T=conv(T1,T2);
T_mul=poly2sym(T)

[A_coe, A_r]=deconv(T1,T3);
T_coe=poly2sym(A_coe)
T_rem=poly2sym(A_r)

********************************************************************************************
page 7

例4.1.1-4 多项式求值,求上式f1(x)在x=0.5处的函数值
T1=[2,5,0,4,1,4];
x=0.5;
y=polyval(T1,x)

*******************************************************************************************

page 8

例4.1.1-5 求多项式f1(x)=2x5+5x4+4x2+x+4的根

T1=[2,5,0,4,1,4];
root=roots(T1)

例4.1.1-6 对多项式f1(x)=2x5+5x4+4x2+x+4进行导数运算

T1=[2,5,0,4,1,4];
h=polyder(T1);
poly2sym(h)

*******************************************************************************************

page 10

例4.1.1-7 对下面给定数据对(x0,y0)求三次拟合多项式，并图示拟合情况。
x0=0:0.1:1;
y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];


x0=0:0.1:1;
y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];
n=3; % 设定拟合次数为3
[p,s]=polyfit(x0,y0,n); % 得到拟合多项式向量和相关偏差信息
f=poly2sym(p);
xx=0:0.01:1;
yy=polyval(p,xx); % 按拟合曲线计算采样值
n1=6; % 设定拟合次数为6
[p1,s1]=polyfit(x0,y0,n1);
yy1=polyval(p1,xx);
plot(xx,yy,'-b',xx,yy1,'-m',x0,y0,'.r','MarkerSize',20); %绘图

% hold on;
y3=polyval(p,0.5);
y6=polyval(p1,0.5);
text(0.5, y3-0.3, 'n=3');
text(0.5, y6+0.2, 'n=6');

///////////////////////////////////////////////////////////////////////////

page 11

例4.1.1-8 按上例所给数据，研究插值，并观察插值和逆合的区别
x0=0:0.1:1;
y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];
xi=0:0.02:1;
yi=interp1(x0,y0,xi,'cubic');
subplot(1,2,1), plot(xi,yi,'-b',x0,y0,'.r','MarkerSize',20);
xlabel('x'), ylabel('y');
[p,s]=polyfit(x0,y0,3);
xx=0:0.01:1;
yy=polyval(p,xx);
subplot(1,2,2);
plot(xx,yy,'-r', x0,y0,'.r','MarkerSize',20);

xlabel('x');
ylabel('y');

****************************************************************************************
page 15
% 直线拟合



format short
format compact
xi=0:0.02:1;
yi1=interp1(x0,y0,xi,'linear');
subplot(1,2,2), plot(xi,yi1,'-b',x0,y0,'.r','MarkerSize',20), xlabel('x')



拟合的使用


format short
format compact
i0=0:0.5:10;
v0=[0,98,202,306,395,506,592,708,789,890,1009,1108,1186,1310,1391,1510,1612,1716,1782,1920,2014];

n1=1;
[p1,s]=polyfit(i0,v0,n1);

n2=2;
[p2,s]=polyfit(i0,v0,n2);

n3=3;
[p3,s]=polyfit(i0,v0,n3);


ii=0:0.1:10;
v1=polyval(p1,ii);
v2=polyval(p2,ii);
v3=polyval(p3,ii);

subplot(1,3,1),plot(ii,v1,'-b',i0,v0,'.r','MarkerSize',8), xlabel('一次拟合');
subplot(1,3,2),plot(ii,v2,'-b',i0,v0,'.r','MarkerSize',8), xlabel('二次拟合');
subplot(1,3,3),plot(ii,v3,'-b',i0,v0,'.r','MarkerSize',8), xlabel('三次拟合');

************************************************************************************************
page 17

A=[2,3,3;4,7,7;-2,4,5];
b=[3,1,-7]';
x=A/b, y=b'/A'

************************************************************************************************
page 12

A=[2,-1,-1;3,4,-2;3,-2,4];
b=[4;11;11];
detA=det(A),r_A=rank(A), r_A_b=rank([A,b])
x=A\b

*************************************************************************************************
page 19

A=[1,2;2,3;3,4];
b=[1;2;3];
x1=inv(A'*A)*A'*b;
x2=A\b

**************************************************************************************************

page 20

A=[1,2,3;2,3,4];
b=[1;2];
x1=A\b;
x2=pinv(A)*b

**************************************************************************************************
page 22

例4.2.3-1 对于超定方程y=Ax, 进行三种解法比较
A=gallery(5);
A(:,1)=[];
y=[1.7, 7.5, 6.3, 0.83, -0.082]';
x=inv(A'*A)*A'*y, xx=pinv(A)*y, xxx=A\y;

nx=norm(x), nxx=norm(xx), nxxx=norm(xxx)

e=norm(y-A*x), ee=norm(y-A*xx), eee=norm(y-A*xxx)

***************************************************************************************************
page 19
A=[2,-1,-1;3,4,-2;3,-2,4]

d=eig(A)
[V,D]=eig(A); A*V V*D
[V,J]=jordan(A)
c=condeig(A)
[V,D,C]=condeig(A)


**********************************************************************************************
page 21

A=[1,2,3;4,5,6;7,8,9]
各列最大值：max(A)
各列最大值：

********************************************************************************************
page 28

例4.3-1 基本统计示例

randn('state',0);
A=randn(1000,4);
AMAX=max(A), AMIN=min(A)
AMED=median(A)
AMEAN=mean(A)
ASTD=std(A)

**********************************************************************************************
rand('state',1)
X=rand(10,3); Y=rand(10,3);
mx=mean(X); Xmx=X-ones(size(X))*diag(mx);
CCX=Xmx'*Xmx/(size(Xmx,1)-1)
CX=cov(X), CY=cov(Y)
Cxy=cov(X,Y)
PX=corrcoef(X)
Pxy=corrcoef(X,Y)


*****************************************************

*******************************************
page 31

dt=0.001;
t=0:dt:2;
y=3*t.^2;
s1=dt*sum(y);
s2=trapz(t,y)
s=dt*cumsum(y);
s3=s(end);

***************************************************************************************************
page 36

(1) 绘制函数图形

x=-10:0.1:10;
y=2.^x-10*x;
plot(x,y)
plot(x,zeros(size(x)),'k');
hold on;
plot(x,y), xlabel('x'), ylabel('y'),title('曲线方程式为：y=2^x-10*x');
hold off;


(2) 获取初步近似解

zoom on;
[tt,yy]=ginput(2);
zoom off;


(3) 利用fzero指令求精确解

func='2^x-10*x';
[x01, err1, exitflag1]=fzero(func, tt(1), []);
[x02, err2, exitflag2]=fzero(func, tt(2), []);

str=['方程的根为: x1= ', num2str(x01), ' 和 x2= ', num2str(x02)];
disp(str);

****************************************************************************************************

page 30

func='2^x-10*x';
x1=0;
x2=8;

[Xmin,Ymin]=fminbnd(func,x1,x2)


保存为func.m
function y=f(x)
y=2.^x-10*x;



x1=0;
x2=8;
Hf=@func;
[X,Ymin]=fminbnd(Hf,x1,x2)

****************************************************************************************************
page 31

例4.4.3-1求f(x)=3x*x在区间[0,2]的积分

y='3*x^2';
a=0;
b=2;
q1=quad(y,a,b)
q2=quadl(y,a,b)



******************************************************************************************************

page 33

例4.5-1 解微分方程：y’=3x2, y(0)=1,微分区间[0,5]

% dy_file.m
function dy=f(x,y)
dy=3*x^2

y0=1;
xspan=[0,5];
[x, y]=ode45('dy_file',xspan,y0);
plot(x,y,'r'), xlabel('x'), ylabel('y')

************************************************
例4.5-2 求微分方程 x''-a.(1-x2)x'+x=0在初始条件x(0)=1,x'(0)=0情况下的解，区间[0,30], a=2


编写M函数文件

% dydt.m
function ydot=dydt(t,y)
a=2;
ydot=[y(2),a*(1-y(1)^2)*y(2)-y(1)]';



解算微分方程

tspan=[0,30];
y0=[1;0];
[tt,yy]=ode45(@dydt,tspan,y0);
plot(tt,zeros(size(tt)),'k'), hold on;
plot(tt,yy(:,1)), title('x(t)'), xlabel('t'), ylabel('x')





************************************************************************************************

page 35

滤波

clear;
randn('state',1);
ws=1000;
t=0:1/ws:0.4;
x=sin(2*pi*10*t)+cos(2*pi*100*t)+0.2*randn(size(t));
wn=ws/2;
[B,A]=butter(10,30/wn);
y=filter(B, A, x);
plot(t,x,'b-',t,y,'r.','MarkerSize',10);