第一篇：單純形法matlab程序

算法實現與分析

算法1．單純形法 具體算例：

minz=?3x1+x2+2x3 3x1+2x2?3x3=6 x1?2x2+x3+x5=4

x1,x2,x3≥0標準化后：

min z=?3x1+x2+2x3+Mx4+Mx5

3x1+2x2?3x3+x4=6 x1?2x2+x3+x5=4

x1,x2,x3,x4,x5≥0用單純形法求解，程序如下： clear clc

M=1000000;

A=[3,2,-3,1,0;1,-2,1,0,1];%系數矩陣 C=[-3,1,2,M,M,0];%價值矩陣 B=[6;4];Xt=[4 5];

for i=1:length(C)-1 D=0;

for j=1:length(Xt)

D=D+A(j,i)*C(Xt(j));

end

xi(i)=C(i)-D;end s=[];

for i=1:length(xi)

if xi(i)<0 s=[s,i];

end end

f=length(s);h=1;

while(f)

for k=1:length(s)j=1;A x=[];

for i=1:length(Xt)

if A(i,s(k))>0 x(j)=i;j=j+1;

end end x

if(length(x)+1==1)break;end y=1 x

for i=1:length(x)

if B(x(i))/A(x(i),s(k))

end end y=x(y);end

y1=Xt(y);%??3?±?á? s k

aa=A(y,s(k))%s(k)?a??è?±?á? A(y,:)=A(y,:)./aa;B(y,:)=B(y,:)./aa;z=[];

for i=1:length(Xt)z=[z,i];end

z z(y)=[];z Xt

for i=1:length(z);yz=-A(z(i),s(k))

A(z(i),:)=A(z(i),:)+A(y,:).*yz B(z(i))B(y)yz

B(z(i))=B(z(i))+B(y).*yz end

for i=1:length(Xt)

if Xt(i)==y1 Xt(i)=s(k);break

end end Xt

disp('×a??oó')A=A B=B AB=[A,B];

for i=1:length(C)D=0;

for j=1:length(Xt)D=D+AB(j,i)*C(Xt(j));

end

xi(i)=C(i)-D;

end xi s=[];

for i=1:length(xi)-1

if xi(i)<0 s=[s,i];

end

end s

vpa([A,B;C]);f=length(s);h=h+1;

if h==5

break

end end

-xi(length(xi))

第二篇：改進單純形法matlab程序

clear clc

X=[1 2 3 4 5];A=[ 1 2 1 0 0;4 0 0 1 0;0 4 0 0 1];C=[2 3 0 0 0 ];b=[8;16;12];t=[3 4 5];B0=A(:,t);while 1

CB0=C(:,t);XN01=X;

for i=1:length(t);

for j=1:length(X);

if XN01(j)==t(i)

XN01(j)=0;

end

end

end j=1;

for i=1:length(X);

if XN01(i)~=0

XN0(j)=XN01(i);

j=j+1;

end

end

for j=1:length(XN0);

CN0(j)=C(XN0(j));

end N0=[];

for i=1:length(XN0);

N0=[N0,A(:,XN0(i))];

end

xiN0=CN0-CB0*B0*N0;j=1;z=[];

for i=1:length(xiN0)

if xiN0(i)>0

z(j)=i;

j=j+1;

end

end

if length(z)+1==1;

break;

end n=1;

for i=1:length(z)

if z(i)>z(n)

n=i;

end

end

k=XN0(z(n));%換入變量 B=B0*b;

P=B0*A(:,k);j=1;

for i=1:length(P)

if P(i)>0

x(j)=i;

j=j+1;

end

end y=1;

for i=1:length(x)

if B(x(y))/P(x(y))>B(x(i))/P(x(i))

y=i;

end

end

y1=x(y);

y=t(y1);%換出變量

for i=1:length(t)

if t(i)==y

m=i;

break;

end

end

t(m)=k;

P2=B0*A(:,k);q=P2(y1);P2(y1)=-1;P2=-P2./q;

E=[1 0 0;0 1 0;0 0 1];E(:,m)=P2;B0=E*B0;end

CB0*B0*b

第三篇：運籌學單純形法matlab程序

function [xx,fm]=myprgmh(m,n,A,b,c)B0=A(:,1:m);cb=c(:,1:m);xx=1:n;sgm=c-cb*B0^-1*A;h=-1;sta=ones(m,1);for i=m+1:n

if sgm(i)>0

h=1;

end end

while h>0

[msg,mk]=max(sgm);

for i=1:m

sta(i)=b(i)/A(i,mk);

end

[mst,mr]=min(sta);

zy=A(mr,mk);

for i=1:m

if i==mr

for j=1:n

A(i,j)=A(i,j)/zy;

end

b(i)=b(i)/zy;

end

end

for i=1:m

if i~=mr

for j=1:n

A(i,j)=A(i,j)-A(i,mk)*A(mr,j);

end

b(i)=b(i)-A(i,mk)*b(mr);

end

end

B1=A(:,1:m);

cb(mr)=c(mk);

xx(mr)=mk;

sgm=c-cb*B1*A;

for i=m+1:n

if sgm(i)>0

h=1;

end

end

end fm=c*xx;

第四篇：線性規劃單純形法matlab解法

線性規劃單純形法matlab解法

%單純形法matlab程序-ssimplex % 求解標準型線性規劃:max c*x;s.t.A*x=b;x>=0 % 本函數中的A是單純初始表，包括:最后一行是初始的檢驗數，最后一列是資源向量b % N是初始的基變量的下標

% 輸出變量sol是最優解, 其中松弛變量（或剩余變量）可能不為0 % 輸出變量val是最優目標值，kk是迭代次數 % 例：max 2*x1+3*x2 % s.t.x1+2*x2<=8 % 4*x1<=16 % 4*x2<=12 % x1,x2>=0 % 加入松馳變量，化為標準型，得到 % A=[1 2 1 0 0 8;% 4 0 0 1 0 16;% 0 4 0 0 1 12;% 2 3 0 0 0 0];% N=[3 4 5];% [sol,val,kk]=ssimplex(A,N)% 然后執行 [sol,val,kk]=ssimplex(A,N)就可以了。function [sol,val,kk]=ssimplex(A,N)[mA,nA]=size(A);kk=0;% 迭代次數 flag=1;while flag kk=kk+1;if A(mA,:)<=0 % 已找到最優解 flag=0;sol=zeros(1,nA-1);for i=1:mA-1 sol(N(i))=A(i,nA);end val=-A(mA,nA);else for i=1:nA-1 if A(mA,i)>0&A(1:mA-1,i)<=0 % 問題有無界解 disp('have infinite solution!');flag=0;break;end end if flag % 還不是最優表，進行轉軸運算 temp=0;for i=1:nA-1 if A(mA,i)>temp temp=A(mA,i);inb=i;% 進基變量的下標 end end sita=zeros(1,mA-1);for i=1:mA-1 if A(i,inb)>0 sita(i)=A(i,nA)/A(i,inb);end end temp=inf;for i=1:mA-1 if sita(i)>0&sita(i)

A(outb,:)=A(outb,:)/A(outb,inb);for i=1:mA if i~=outb A(i,:)=A(i,:)-A(outb,:)*A(i,inb);End End End End end

第五篇：MATLAB程序總結

%******************************讀取數據************************************* MATLAB讀取數據

Xlsread(‘lujing’,’mingcheng’)%******************************輸出數據************************************* MATLAB輸出Excel數據

Xlswrite(‘eular.xla’,[A],’sheet 1’)%******************************三維曲面************************************* 三維曲面

x1=[0.05 0.05 0.05 0.1 0.05 0.1 0.2 0 0 0 0 0 0 0 0];x2=[0.1 0.25 0.12 0.12 0 0 0 0 0.05 0.08 0.1 0.12 0.15 0.2 0.25];y1=[153 130 146 133 160 154 140 133 165 169 171 186 175 156 152];y2=[7.22 5.57 6.66 6 6.6 8.21 5.1 4.66 7.58 8.26 8.99 8.73 7.71 7 6.49];y3=[112 100 131 88 80 72 72 80 116 83 80 80 64 68 74];scatter3(x1,x2,y2)figure [X,Y,Z]=griddata(x1,x2,y2,linspace(min(x1),max(x1))',linspace(min(x2),max(x2)),'v4');%插值 pcolor(X,Y,Z);shading interp%偽彩色圖 figure,contourf(X,Y,Z)%等高線圖 figure,surf(X,Y,Z);%三維曲面

%******************************各種距離************************************* pdist(X)— 樣本X中各n維向量的歐氏距離 pdist(X, 'cityblock')— 各n維向量的絕對距離 pdist(X, ‘minkowski'，r)— 閔可夫斯基距離 pdist(X, 'mahal')— 各n維向量的馬氏距離 mean(A)— A中各列向量的均值

Var(A)— A中各列向量的方差

Std(A)— A中各列向量的標準差

Cov(A)— A中各列向量的協方差矩陣 Corrcoef(A)— A中各列向量的相關矩陣 命令：[x,fval] = fminbnd(F,a,b)%******************************泰勒公式************************************* 計算F在a,b之間取極小值時的x與y(即fval).格式：taylor(F,a)功能：F在x=a處的泰勒級數前5項 格式：taylor(F,v)功能：F對變量v的泰勒展式前5項

格式：taylor(F,v,n)功能：求F的n 階泰勒展式,且

(n缺省時默認為 5）

%******************************一維插值************************************* x=[-2*pi:0.5*pi:2*pi];y=cos(x);xi=-2*pi:0.3*pi:2*pi;y_nearest=interp1(x,y, xi,‘nearest’)%最近鄰插值 y_linear= interp1(x,y,xi)%雙線性差值

y_spline= interp1(x,y,xi, ‘spline’)%三次樣條插值 y_cubic= interp1(x,y,xi, ‘cubic’)%雙三次插值

plot(x,y,'o',xi,y_nearest,'-',xi,y_linear, 'r* ', xi,y_spline,'k:',xi,y_cubic,'k-');legend('original data','nearest','linear','spline','cubic')%******************************曲線擬合************************************* x1=[2:16];y1=[6.42,8.2,9.58,9.5,9.7,10,9.93,9.99,10.49,10.59,10.6,10.8,10.6,10.9,10.76];b01=[0.1435,0.084];

%初始參數值

fun1=inline('x./(b(1)+b(2)*x)','b','x');% 定義函數 [b1,r1,j1]=nlinfit(x1,y1,fun1,b01);y=x1./(0.1152+0.0845*x1);%根據b1寫出具體函數

plot(x1,y1,'*',x1,y,'-or');R2=1-sum((y-y1).^2)/sum((y1-mean(y1)).^2)%*****************************殘差與置信區間******************************** X=[ones(30,1), x1‘,x2’,x3‘];

%輸入自變量（注意1與自變量組成的矩陣）[b,bint,r,rint,s]=regress(y',X);%多元線性回歸 b,bint,s rcoplot(r,rint)%殘差及其置信區間作圖

a=[1,3:9,11:30];X1=X(a,:);y1=y(a);%剔除異常點

[b1,bint1,r1,rint1,s1]=regress(y1',X1);b1,bint1,s1 %再次線性回歸 rcoplot(r1,rint1)

%------------%

Copula理論及其在matlab中的實現程序應用實例 %------------

%******************************讀取數據************************************* % 從文件hushi.xls中讀取數據 hushi = xlsread('hushi.xls');% 提取矩陣hushi的第5列數據，即滬市的日收益率數據 X = hushi(:,5);% 從文件shenshi.xls中讀取數據 shenshi = xlsread('shenshi.xls');% 提取矩陣shenshi的第5列數據，即深市的日收益率數據 Y = shenshi(:,5);

%****************************繪制頻率直方圖*********************************

% 調用ecdf函數和ecdfhist函數繪制滬、深兩市日收益率的頻率直方圖 [fx, xc] = ecdf(X);figure;ecdfhist(fx, xc, 30);xlabel('滬市日收益率');% 為X軸加標簽 ylabel('f(x)');% 為Y軸加標簽 [fy, yc] = ecdf(Y);figure;ecdfhist(fy, yc, 30);xlabel('深市日收益率');% 為X軸加標簽 ylabel('f(y)');% 為Y軸加標簽

%****************************計算偏度和峰度********************************* % 計算X和Y的偏度 xs = skewness(X)ys = skewness(Y)

% 計算X和Y的峰度 kx = kurtosis(X)ky = kurtosis(Y)

%******************************正態性檢驗*********************************** % 分別調用jbtest、kstest和lillietest函數對X進行正態性檢驗 [h,p] = jbtest(X)% Jarque-Bera檢驗

[h,p] = kstest(X,[X,normcdf(X,mean(X),std(X))])% Kolmogorov-Smirnov檢驗 [h, p] = lillietest(X)% Lilliefors檢驗

% 分別調用jbtest、kstest和lillietest函數對Y進行正態性檢驗 [h,p] = jbtest(Y)% Jarque-Bera檢驗

[h,p] = kstest(Y,[Y,normcdf(Y,mean(Y),std(Y))])% Kolmogorov-Smirnov檢驗 [h, p] = lillietest(Y)% Lilliefors檢驗

%****************************求經驗分布函數值******************************* % 調用ecdf函數求X和Y的經驗分布函數 [fx, Xsort] = ecdf(X);[fy, Ysort] = ecdf(Y);% 調用spline函數，利用樣條插值法求原始樣本點處的經驗分布函數值 U1 = spline(Xsort(2:end),fx(2:end),X);V1 = spline(Ysort(2:end),fy(2:end),Y);

% 調用ecdf函數求X和Y的經驗分布函數 [fx, Xsort] = ecdf(X);[fy, Ysort] = ecdf(Y);% 提取fx和fy的第2個至最后一個元素，即排序后樣本點處的經驗分布函數值 fx = fx(2:end);fy = fy(2:end);

% 通過排序和反排序恢復原始樣本點處的經驗分布函數值U1和V1 [Xsort,id] = sort(X);[idsort,id] = sort(id);U1 = fx(id);[Ysort,id] = sort(Y);[idsort,id] = sort(id);V1 = fy(id);

%*******************************核分布估計********************************** % 調用ksdensity函數分別計算原始樣本X和Y處的核分布估計值 U2 = ksdensity(X,X,'function','cdf');V2 = ksdensity(Y,Y,'function','cdf');

% **********************繪制經驗分布函數圖和核分布估計圖********************** [Xsort,id] = sort(X);% 為了作圖的需要，對X進行排序 figure;% 新建一個圖形窗口

plot(Xsort,U1(id),'c','LineWidth',5);% 繪制滬市日收益率的經驗分布函數圖 hold on plot(Xsort,U2(id),'k-.','LineWidth',2);% 繪制滬市日收益率的核分布估計圖 legend('經驗分布函數','核分布估計', 'Location','NorthWest');% 加標注框 xlabel('滬市日收益率');% 為X軸加標簽 ylabel('F(x)');% 為Y軸加標簽

[Ysort,id] = sort(Y);% 為了作圖的需要，對Y進行排序 figure;% 新建一個圖形窗口

plot(Ysort,V1(id),'c','LineWidth',5);% 繪制深市日收益率的經驗分布函數圖 hold on plot(Ysort,V2(id),'k-.','LineWidth',2);% 繪制深市日收益率的核分布估計圖 legend('經驗分布函數','核分布估計', 'Location','NorthWest');% 加標注框 xlabel('深市日收益率');% 為X軸加標簽 ylabel('F(x)');% 為Y軸加標簽

%****************************繪制二元頻數直方圖***************************** % 調用ksdensity函數分別計算原始樣本X和Y處的核分布估計值 U = ksdensity(X,X,'function','cdf');V = ksdensity(Y,Y,'function','cdf');figure;% 新建一個圖形窗口

% 繪制邊緣分布的二元頻數直方圖，hist3([U(:)V(:)],[30,30])xlabel('U（滬市）');% 為X軸加標簽 ylabel('V（深市）');% 為Y軸加標簽 zlabel('頻數');% 為z軸加標簽

%****************************繪制二元頻率直方圖***************************** figure;% 新建一個圖形窗口

% 繪制邊緣分布的二元頻數直方圖，hist3([U(:)V(:)],[30,30])h = get(gca, 'Children');% 獲取頻數直方圖的句柄值 cuv = get(h, 'ZData');% 獲取頻數直方圖的Z軸坐標

set(h,'ZData',cuv*30*30/length(X));% 對頻數直方圖的Z軸坐標作變換 xlabel('U（滬市）');% 為X軸加標簽 ylabel('V（深市）');% 為Y軸加標簽 zlabel('c(u,v)');% 為z軸加標簽

%***********************求Copula中參數的估計值****************************** % 調用copulafit函數估計二元正態Copula中的線性相關參數 rho_norm = copulafit('Gaussian',[U(:), V(:)])% 調用copulafit函數估計二元t-Copula中的線性相關參數和自由度 [rho_t,nuhat,nuci] = copulafit('t',[U(:), V(:)])

%********************繪制Copula的密度函數和分布函數圖************************ [Udata,Vdata] = meshgrid(linspace(0,1,31));% 為繪圖需要，產生新的網格數據 % 調用copulapdf函數計算網格點上的二元正態Copula密度函數值 Cpdf_norm = copulapdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);% 調用copulacdf函數計算網格點上的二元正態Copula分布函數值 Ccdf_norm = copulacdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);% 調用copulapdf函數計算網格點上的二元t-Copula密度函數值 Cpdf_t = copulapdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);% 調用copulacdf函數計算網格點上的二元t-Copula分布函數值 Ccdf_t = copulacdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);% 繪制二元正態Copula的密度函數和分布函數圖 figure;% 新建圖形窗口

surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata)));% 繪制二元正態Copula密度函數圖 xlabel('U');% 為X軸加標簽 ylabel('V');% 為Y軸加標簽 zlabel('c(u,v)');% 為z軸加標簽 figure;% 新建圖形窗口

surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata)));% 繪制二元正態Copula分布函數圖 xlabel('U');% 為X軸加標簽 ylabel('V');% 為Y軸加標簽 zlabel('C(u,v)');% 為z軸加標簽

% 繪制二元t-Copula的密度函數和分布函數圖 figure;% 新建圖形窗口

surf(Udata,Vdata,reshape(Cpdf_t,size(Udata)));% 繪制二元t-Copula密度函數圖 xlabel('U');% 為X軸加標簽 ylabel('V');% 為Y軸加標簽 zlabel('c(u,v)');% 為z軸加標簽 figure;% 新建圖形窗口

surf(Udata,Vdata,reshape(Ccdf_t,size(Udata)));% 繪制二元t-Copula分布函數圖 xlabel('U');% 為X軸加標簽 ylabel('V');% 為Y軸加標簽 zlabel('C(u,v)');% 為z軸加標簽

%**************求Kendall秩相關系數和Spearman秩相關系數*********************** % 調用copulastat函數求二元正態Copula對應的Kendall秩相關系數 Kendall_norm = copulastat('Gaussian',rho_norm)% 調用copulastat函數求二元正態Copula對應的Spearman秩相關系數 Spearman_norm = copulastat('Gaussian',rho_norm,'type','Spearman')% 調用copulastat函數求二元t-Copula對應的Kendall秩相關系數 Kendall_t = copulastat('t',rho_t)% 調用copulastat函數求二元t-Copula對應的Spearman秩相關系數 Spearman_t = copulastat('t',rho_t,'type','Spearman')

% 直接根據滬、深兩市日收益率的原始觀測數據，調用corr函數求Kendall秩相關系數 Kendall = corr([X,Y],'type','Kendall')% 直接根據滬、深兩市日收益率的原始觀測數據，調用corr函數求Spearman秩相關系數 Spearman = corr([X,Y],'type','Spearman')

%******************************模型評價************************************* % 調用ecdf函數求X和Y的經驗分布函數 [fx, Xsort] = ecdf(X);[fy, Ysort] = ecdf(Y);% 調用spline函數，利用樣條插值法求原始樣本點處的經驗分布函數值 U = spline(Xsort(2:end),fx(2:end),X);V = spline(Ysort(2:end),fy(2:end),Y);% 定義經驗Copula函數C(u,v)C = @(u,v)mean((U <= u).*(V <= v));% 為作圖的需要，產生新的網格數據 [Udata,Vdata] = meshgrid(linspace(0,1,31));% 通過循環計算經驗Copula函數在新產生的網格點處的函數值 for i=1:numel(Udata)

CopulaEmpirical(i)= C(Udata(i),Vdata(i));end

figure;% 新建圖形窗口

% 繪制經驗Copula分布函數圖像

surf(Udata,Vdata,reshape(CopulaEmpirical,size(Udata)))xlabel('U');% 為X軸加標簽 ylabel('V');% 為Y軸加標簽

zlabel('Empirical Copula C(u,v)');% 為z軸加標簽

% 通過循環計算經驗Copula函數在原始樣本點處的函數值 CUV = zeros(size(U(:)));for i=1:numel(U)

CUV(i)= C(U(i),V(i));end

% 計算線性相關參數為0.9264的二元正態Copula函數在原始樣本點處的函數值 rho_norm = 0.9264;Cgau = copulacdf('Gaussian',[U(:), V(:)],rho_norm);% 計算線性相關參數為0.9325，自由度為4的二元t-Copula函數在原始樣本點處的函數值 rho_t = 0.9325;k = 4.0089;Ct = copulacdf('t',[U(:), V(:)],rho_t,k);% 計算平方歐氏距離

dgau2 =(CUV-Cgau)'*(CUV-Cgau)dt2 =(CUV-Ct)'*(CUV-Ct)

%******************************灰色預測*************************************

A=[ 96 144 194 276 383 466 554 652 747 832 972

169 235 400 459 565 695 805 881 1011 1139

151 238 335 425 541 641 739 866 975 1087 1238

164 263 376 531 600 711 913 1038 1173 1296 1497

182 318 445 576 708 856 1000 1145 1292 1435 1667

216 361 504 642 818 979 1142 1305 1479 1644 1920] m=size(A,2);x0=mean(A,2);x1=cumsum(x0);alpha=0.4;n=length(x0);z1=alpha*x1(2:n)+(1-alpha)*x1(1:n-1);Y=x0(2:n);B=[-z1,ones(n-1,1)];ab=BY

x_A=(x0(1)-ab(2)/ab(1))*(exp(-ab(1)*n)-exp(-ab(1)*(n-1)))

圖標參數************************************* y 黃-實線.點< 小于號 m 紫： 點線o 圓s 正方形 c 青-.點劃線x 叉號d 菱形 r 紅--虛線+ 加號h 六角星 g 綠* 星號p 五角星 b 藍v 向下三角形 w 白^ 向上三角形 k 黑> 大于號

%******************************圖形修飾*************************************

函數 含義

grid on(/off)給當前圖形標記添加（取消）網絡 xlable(‘string’)標記橫坐標 ylabel(‘string’)標記縱坐標 title(‘string’)給圖形添加標題

text(x,y,’string’)在圖形的任意位置增加說明性文本信息 gtext(‘string’)利用鼠標添加說明性文本信息

axis([xmin xmax ymin ymax])設置坐標軸的最小最大值

例5.1.2 給例5.1.1 的圖形中加入網絡和標記。（見圖5.1.3 和5.1.4）>> x=0:pi/10:2*pi;>> y1=sin(x);>> y2=cos(x);>> plot(x,y1,x,y2)>> grid on >> xlabel('independent variable X')>> ylabel('Dependent Variable Y1 & Y2')>> title('Sine and Cosine Curve')>> text(1.5,0.3,'cos(x)')>> gtext('sin(x)')>> axis([0 2*pi-0.9 0.9])

圖5.1.3 使用了圖形修飾的plot 函數繪制的正弦曲線

例5.2.1 繪制方程x=t y=sin(t)z=cos(t)

在t=[0,2*pi]上的空間方程。（見圖5.2.1）

>> clf >> x=0:pi/10:2*pi;>> y1=sin(x);>> y2=cos(x);>> plot3(y1,y2,x,'m:p')>> grid on >> xlabel('Dependent Variable Y1')>> ylabel('Dependent Variable Y2')>> zlabel('Independent Variable X')>> title('Sine and Cosine Curve')

圖5.2.1 函數plot 繪制的三維曲線圖

%******************************三維圖形*************************************

例5.2.3 繪制方程

sin((x^2+y^2)^(1/2))z =---------------------(x^2+y^2)^(1/2)

在x∈[-7.5,7.5];y∈[-7.5,7.5] 的圖形。

>> x=-7.5:0.5:7.5;y=x;>> [X,Y]=meshgrid(x,y);>> R=sqrt(X.^2+Y.^2)+eps;>> Z=sin(R)./R;>> surf(X,Y,Z)>> xlabel('X 軸方向')>> ylabel('Y 軸方向')>> zlabel('Z 軸方向')(見圖5.2.4)

圖5.2.4

例5.2.4 繪制由方程形成的立體圖。（見圖5.2.5）z=x*exp(-(x^2+y^2))

>> clear >> x=-2:0.1:2;y=x;>> [X,Y]=meshgrid(x,y);>> Z=X.*exp(-X.^2-Y.^2);>> surf(X,Y,Z)

圖5.2.5

%******************************三維多角度觀察*********************************

MTALAB 允許用戶設置觀察點，其指令是： view(azimuth,elevation)其中方位角azimuth 是觀察點和坐標原點連線在x-y平面的投影和y 軸負方向的夾角，仰

角elevation 是觀察點與坐標原點的連線和x-y平面的夾角。對于這兩個角度，三維圖形的

默認值分別是-37.5 和30，二維圖形的默認值是0 和90。

例5.2.5 從不同的角度觀察高斯矩陣的曲面。

>> z=peaks(40);>> subplot(2,2,1);>> mesh(z);>> subplot(2,2,2);>> mesh(z);>> view(-37.5,-30);>> subplot(2,2,3);>> mesh(z);>> view(180,0);>> subplot(2,2,4);>> mesh(z);>> view(0,90);

圖5.2.6 對應不同觀察點的三維曲面圖

%******************************其他圖形*************************************

除了plot 繪圖函數以外，在有些場合對繪制的曲線會有一些特殊要求，這就要其他函

數來實現，常用的幾種函數如下（見表5.3.1）

表5.3.1 其他圖形函數表 函數 含義

loglog

使用對數坐標系繪圖

semilogx

橫坐標為對數坐標軸，縱坐標為線性坐標軸 semilogy

橫坐標為線性坐標軸，縱坐標為對數坐標軸 polar

繪制極坐標圖 fill

繪制實心圖 bar

繪制直方圖 pie

繪制餅圖 area

繪制面積圖 quiver 繪制向量場圖 stairs 繪制階梯圖 sterm 繪制火柴桿圖

例5.3.1 >> x=0:pi/10:2*pi;>> y1=sin(x);>> subplot(2,2,1);>> plot(x,y1);>> subplot(2,2,2);>> bar(x,y1);>> subplot(2,2,3);>> fill(x,y1,'g');>> subplot(2,2,4);>> stairs(x,y1,'k');

%******************************直方圖*************************************

函數bar(x)可以繪制直方圖，這對統計或者數據采集非常直觀實用。它共有四種形式： bar,bar3,barh 和bar3h，其中bar 和bar3 分別用來繪制二維和三維豎直方圖，barh 和b ar3h 分別用來繪制二維和三維水平直方圖，調用格式是：

bar(x,y)其中x 必須單調遞增或遞減，y 為n m× 矩陣，可視化結果為m 組，每

組n 個垂直柱，也就是把y 的行畫在一起，同一列的數據用相同的顏色表示； bar(x,y,width)（或bar(y,width)）指定每個直方條的寬度，如width>1，則直方條會重

疊，默認值為width=0.8；

bar(…,’grouped’)使同一組直方條緊緊靠在一起； bar(…,’stack’)把同一組數據描述在一個直方條上。

例5.3.2

>> y=[5 3 2 9;4 7 2 7;1 5 7 3];>> subplot(2,2,1),bar(y)>> x=[5 9 11];>> subplot(2,2,2),bar3(x,y)>> subplot(2,2,3),bar(x,y,'grouped')>> subplot(2,2,4),bar(rand(2,3),.75,'stack')

%******************************面積圖*************************************

5.3.2 面積圖

函數area 用來繪制面積圖，面積圖在plot 的基礎上填充x 軸和曲線之間的面積，該圖

用于查看某個數在該列所有數的總和中所占的比例。

例5.3.3

>> x=-3:3;>> y=[3 2 5;6 1 8;7 4 9;6 3 7;8 2 9;4 2 9;3 1 7];>> area(x,y)

%******************************餅圖*************************************

5.3.3 餅圖

函數pie 用來繪制餅圖，它可以形象地表示出向量中各元素所占比例。其調用格式是：

pie(x)x 中的元素通過x/sum(x)進行歸一化，以確定餅圖中的份額； pie(x,explode)向量explode 和x 元素數相同，用來指出需要分開的餅片，explode 中

不為零的部分會被分開。

例5.3.4 設某班的某課程的考試成績如下：90 分以上有32 人，81 至90 有58 人，71 至80 分有27 人，60 至70 分為21 人，60 分以下有16 人，畫出餅圖。（見圖5.3.4）

>> x=[32 58 27 21 16];>> explode0=[1 0 0 0 0];>> subplot(1,2,1)>> pie(x,explode0)>> explode1=[0 0 0 0 1];>> subplot(1,2,2)>> pie(x,explode1)

5.3.4 不同坐標系中的繪圖

Semilogx,semilogy,loglo,polar(theta,rho)的使用方法和plot 完全類似，不同的只是繪

制到不同的圖形坐標上。函數semilogx 繪制x 軸為對數標度的圖形，在半對數坐標系中繪圖； 函數semilogy 繪制y 軸為對數標度的圖形；函數loglog 繪制兩個軸都為對數間隔的圖形；

函數polar(theta,rho)繪制極坐標圖形，其中theta 為相角，rho 為其對應的半徑。

例5.3.5 繪制ρ=acos(3θ),a=2 的圖形。（見圖5.3.5）

>> theta=-pi:pi/80:pi;>> polar(theta,2*cos(3*theta))

圖5.3.5 極坐標圖

5.4 符號表達式繪圖

MATLAB 軟件提供了將表達式進行圖形顯示的功能。完成此功能需調用fplot 函數和

ezplot 函數。

函數fplot 用來繪制數學函數，其調用格式為： fplot(fun,lims)其中fun 就是所要繪制的函數，可以是定義函數的M 文件名，也可以是以x 為變量的可計

算字符串。例如’diric(x,10)’或’[sin(x),cos(x)]’，對于向量x 的每個元素，函數 fun(x)必須返回一個行向量。如果fun 返回[f1(x),f2(x),f3(x)]，輸入[x1;x2]，就會返回矩陣

f1(x1)f2(x1)f3(x1)f1(x2)f2(x2)f3(x2)lims=[XMIN XMAX YMIN YMAX]限定了x,y 軸上的繪圖空間。

例5.4.1 >> subplot(2,2,1),fplot('humps',[0 1])>> subplot(2,2,2),fplot('abs(exp(-j*x*(0:9))*ones(10,1))',[0 2*pi])>> subplot(2,2,3),fplot('[tan(x),sin(x),cos(x)]',2*pi*[-1 1-1 1])>> subplot(2,2,4),fplot('sin(1./x)',[0.01 0.1],1e-3)

圖5.4.1 fplot 函數繪制表達式圖形

ezplot 函數是簡捷繪圖指令之一，它無需數據準備，直接畫出函數圖形，基本調用格式

為ezplot(f)其中f 是字符串或代表數學函數的符號表達式，只有一個符號變量，可以是x，缺省情況下

x 軸的繪圖區域為[-π, π ]，但我們可以用ezplot(f,xmin,xmax)或ezplot(f,[xmin,xmax])來指定x 的范圍。

例5.4.2

>> y='x^2';>> subplot(1,2,1)>> ezplot(y)>> subplot(1,2,2)>> y='sin(x)';>> ezplot(y,[-pi,pi])

圖5.4.2 ezplot 函數繪制表達式圖形

5.5 plot 函數

MATLAB 對數據是按列存儲和計算的，運用plot(x)時，當x 為一個向量時，以其元

素為縱坐標，其序號為橫坐標值繪制曲線。當x 為實矩陣時，則以其序號為橫坐標，按列

繪制每列元素相對于序號的曲線，當x 為n m× 矩陣時，就有n 條曲線。

如果x,y 是同維向量，plot(x,y)指令以x 元素為橫坐標值，y 元素為縱坐標值繪制曲線。

如x 是向量，y 是有一維與x 元素數量相等的矩陣，則以x 為共同橫坐標，按列繪制y 每

列元素值，曲線數為y 的另一維的元素數。如果x,y 是同維矩陣，則以x,y 對應列元素為、縱坐標分別繪制曲線，數目等于矩陣的列數。

例5.5.1

>> x=[3 5 10 8];>> subplot(2,2,1)>> plot(x)>> x=[3 5 10 8;7 2 9 4;2 7 2 7]';>> subplot(2,2,2)>> plot(x)>> x=[3 5 6 8];>> y=[1 5 10 4];>> subplot(2,2,3)>> plot(x,y)>> x=[1 3 5 7;2 4 6 8]';>> y=[6 2 5 10;3 5 2 6]';>> subplot(2,2,4)>> plot(x,y,'k:*')