一维等离子体FDTD的Matlab源代码(两种方法)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% 1D %%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%初始化

clear; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%系统参数

TimeT=3000;%迭代次数

KE=2000;%网格树木

kc=450;%源的位置

kpstart=500;%等离子体开始位置

kpstop=1000;%等离子体终止位置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%物理参数

c0=3e8;%真空中波速

zdelta=1e-9;%网格大小

dt=zdelta/(2*c0);%时间间隔

f=900e12;%Gause脉冲的载频

d=3e-15%脉冲底座宽度

t0=2.25/f;%脉冲中心时间

u0=57e12%碰撞频率

fpe=2000e12;%等离子体频率

wpe=2*pi*fpe;%等离子体圆频率

epsz=1/(4*pi*9*10^9); % 真空介电常数

mu=1/(c0^2*epsz);%磁常数

ex_low_m1=0;

ex_low_m2=0;

ex_high_m1=0;

ex_high_m2=0;

a0=2*u0/dt+(2/dt)^2;

a1=-8/(dt)^2;

a2=-2*u0/dt+(2/dt)^2;

b0=wpe^2+2*u0/dt+(2/dt)^2;

b1=2*wpe^2-8/(dt)^2;

b2=wpe^2-2*u0/dt+(2/dt)^2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%初始化电磁场

Ex=zeros(1,KE);

Ex_Pre=zeros(1,KE);

Hy=zeros(1,KE);

Hy_Pre=zeros(1,KE);

Dx=zeros(1,KE);

Dx_Pre=zeros(1,KE);

Sx1=zeros(1,KE);

Sx2=zeros(1,KE);

Sx3=zeros(1,KE);

Sx=zeros(1,KE);

Dx=Ex;

Dx1=Ex;

Dx2=Ex;

Dx3=Ex;

Ex1=Ex;

Ex2=Ex; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%开始计算

for T=1:TimeT

%%%保存前一时间的电磁场

Ex_Pre=Ex;

Hy_Pre=Hy;

%%%%中间差分计算Dx

for i=2:KE

Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));

end

%%%%%%%%加入源

Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);

%%%计算电场Ex

for i=1:kpstart-1

Ex(i)=Dx(i)/epsz;

end

for i=kpstop+1:KE

Ex(i)=Dx(i)/epsz;

end

Dx3=Dx2;

Dx2=Dx1;

Dx1=Dx;

for i=kpstart:kpstop

Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex 2(i));

end

Ex2=Ex1;

Ex1=Ex;

Sx3=Sx2;

Sx2=Sx1;

Ex(1)=ex_low_m2;

ex_low_m2=ex_low_m1;

ex_low_m1=Ex(2);

Ex(KE)=ex_high_m2;

ex_high_m2=ex_high_m1;

ex_high_m1=Ex(KE-1);

%%%%%%%%%%%%%%%%%%计算磁场

for i=1:KE-1

Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));

end

plot(Ex);

grid on;

pause(0.01);

end

% FDTD Main Function Jobs to Workers %

%*********************************************************************** % 3-D FDTD code with PEC boundaries

%*********************************************************************** %

% This MATLAB M-file implements the finite-difference time-domain

% solution of Maxwell's curl equations over a three-dimensional

% Cartesian space lattice comprised of uniform cubic grid cells.

%

% To illustrate the algorithm, an air-filled rectangular cavity

% resonator is modeled. The length, width, and height of the

% cavity are X cm (x-direction), Y cm (y-direction), and

% Z cm (z-direction), respectively.

%

% The computational domain is truncated using PEC boundary

% conditions:

% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes

% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes

% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes

% These PEC boundaries form the outer lossless walls of the cavity.

%

% The cavity is excited by an additive current source oriented

% along the z-direction. The source waveform is a differentiated

% Gaussian pulse given by

% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),

% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-

% content pulse is approximately 7 GHz. The grid resolution

% (dx = 2 mm) was chosen to provide at least 10 samples per

% wavelength up through 15 GHz.

%

% To execute this M-file, type "fdtd3D" at the MATLAB prompt.

% This M-file displays the FDTD-computed Ez fields at every other

% time step, and records those frames in a movie matrix, M, which

% is played at the end of the simulation using the "movie" command.

%

%*********************************************************************** function [Ex,Ey,Ez]=FDTD3D_Main(handles)

global SimRunStop

% if ~isdir('C:\MATLAB7\work\cavity\figures')

% mkdir 'C:\MATLAB7\work\cavity\figures'

% end

%*********************************************************************** % Grid Partition

%***********************************************************************

p.ip = get(handles.XdirPar,'Value');

p.jp = get(handles.YdirPar,'Value');

p.PN = get(handles.PartNo,'Value');

%*********************************************************************** % Grid Dimensons

%***********************************************************************

ie = get(handles.xslider,'Value'); %number of grid cells in x-direction

je = get(handles.yslider,'Value'); %number of grid cells in y-direction

ke = get(handles.zslider,'Value'); %number of grid cells in z-direction

ib=ie+1;

jb=je+1;

kb=ke+1;

%*********************************************************************** % All Domains Fields Ini.

%***********************************************************************

Ex=zeros(ie,jb,kb);

Ey=zeros(ib,je,kb);

Ez=zeros(ib,jb,ke);

Hx=zeros(ib,je,ke);

Hy=zeros(ie,jb,ke);

Hz=zeros(ie,je,kb);

%*********************************************************************** % Fundamental constants

%*********************************************************************** https://www.wendangku.net/doc/e39209758.html,=2.99792458e8; %speed of light in free space

param.muz=4.0*pi*1.0e-7; %permeability of free space

param.epsz = 1.0/(https://www.wendangku.net/doc/e39209758.html,*https://www.wendangku.net/doc/e39209758.html,*param.muz); %permittivity of free space

%*********************************************************************** % Grid parameters

%***********************************************************************

param.is = get(handles.xsource,'Value'); %location of z-directed current source param.js = get(handles.ysource,'Value'); %location of z-directed current source

param.kobs = floor(ke/2); %Surface of observation

param.dx = get(handles.CellSize,'Value');; %space increment of cubic lattice param.dt = param.dx/(2.0*https://www.wendangku.net/doc/e39209758.html,); %time step

param.nmax = get(handles.TimeStep,'Value');; %total number of time steps

%***********************************************************************

% Differentiated Gaussian pulse excitation

%***********************************************************************

param.rtau=get(handles.tau,'Value')*100e-12;

param.tau=param.rtau/param.dt;

param.ndelay=3*param.tau;

param.Amp = get(handles.sourceamp,'Value')*10e11;

param.srcconst=-param.dt*param.Amp;

%***********************************************************************

% Material parameters

%***********************************************************************

param.eps= get(handles.epsilon,'Value');

param.sig= get(handles.sigma,'Value');

%***********************************************************************

% Updating coefficients

%***********************************************************************

param.ca=(1.0-(param.dt*param.sig)/(2.0*param.epsz*param.eps))/(1.0+(param.dt*param.sig)/( 2.0*param.epsz*param.eps));

param.cb=(param.dt/param.epsz/param.eps/param.dx)/(1.0+(param.dt*param.sig)/(2.0*param.e psz*param.eps));

param.da=1.0;

param.db=param.dt/param.muz/param.dx;

%***********************************************************************

% Calling FDTD Algorithm

%***********************************************************************

ex=zeros(ib,jb,kb);

ey=zeros(ib,jb,kb);

ez=zeros(ib,jb,kb);

hx=zeros(ib,jb,kb);

hy=zeros(ib,jb,kb);

hz=zeros(ib,jb,kb);

[X,Y,Z] = meshgrid(1:ib,1:jb,1:kb); % Grid coordinates

Psim = zeros(param.nmax,1);

Panl = zeros(param.nmax,1);

if ((p.ip == 1)&(p.jp == 0))

x = ceil(ie/p.PN)

param.a = [1:x-1:ie-x];

param.b = [x+1:x-1:ie];

param.c = [1,1];

param.d = [je,je];

m2 = 1;

for n=1:1:param.nmax

for m1=1:1:p.PN

[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);

[hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);

end

[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);

field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);

Dyn_FFT

pause(0.01);

end

elseif ((p.ip == 0)&(p.jp == 1))

y = ceil(je/p.PN);

param.c = [1:y-1:je-y];

param.d = [y+1:y-1:je];

param.a = [1,1];

param.b = [ie,ie];

m1 = 1;

for n=1:1:param.nmax

for m2=1:1:p.PN

[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);

[hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);

end

[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);

field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);

pause(0.01);

end

elseif ((p.ip == 1)&(p.jp == 1))

x = ceil(ie/p.PN);

param.a = [1:x-1:ie-x];

param.b = [x+1:x-1:ie];

y = ceil(je/p.PN);

param.c = [1:y-1:je-y];

param.d = [y+1:y-1:je];

for n=1:1:param.nmax

for m2=1:1:p.PN

for m1=1:1:p.PN

[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);

[hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);

end

end

[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);

field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);

pause(0.01);

end

else

param.a = 1;

param.b = ie;

param.c = 1;

param.d = je;

m1 = 1;m2=1;

for n=1:1:param.nmax

[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);

[hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);

SimRunStop = get(handles.Stop_sim,'value');

if SimRunStop == 1

h = warndlg('Simulation Run is Stopped by User !','!! Warning !!');

waitfor(h);

break;

end

[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);

if n>=2

Panl(n)=Panl(n) + Panl(n-1);

end

field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);

pause(0.01);

end

end

%文件2

% Cavity Field Viz. %

function [] = field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p)

%*********************************************************************** % Cross-Section initialization

%***********************************************************************

tview = squeeze(ez(:,:,param.kobs));

sview = squeeze(ez(:,param.js,:));

ax1 = handles.axes1;

ax2 = handles.axes2;

ax3 = handles.axes3;

%*********************************************************************** % Visualize fields

%***********************************************************************

timestep=int2str(n);

ezview=squeeze(ez(:,:,param.kobs));

exview=squeeze(ex(:,:,param.kobs));

eyview=squeeze(ey(:,:,param.kobs));

xmin = min(X(:));

xmax = max(X(:));

ymin = min(Y(:));

ymax = max(Y(:));

zmin = min(Z(:));

daspect([2,2,1])

xrange = linspace(xmin,xmax,8);

yrange = linspace(ymin,ymax,8);

zrange = 3:4:15;

[cx cy cz] = meshgrid(xrange,yrange,zrange);

% sview=squeeze(ez(:,param.js,:));

rect = [-50 -35 360 210];

rectp = [-50 -40 350 260];

axes(ax1);

axis tight

set(gca,'nextplot','replacechildren');

E_total = sqrt(ex.^2 + ey.^2 + ez.^2);

etview=squeeze(E_total(:,:,param.kobs));

sview = squeeze(E_total(:,param.js,:));

surf(etview');

shading interp;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Total E-Field, time step = ',timestep],'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD');

xlabel('i coordinate');

ylabel('j coordinate');

set(gca,'fontname','Times New Roman','fontsize',10);

% F1 = getframe(gca,rect);

% M1 = frame2im(F1);

% filename = fullfile('C:\MATLAB7\work\cavity\figures',strcat('XY_ETotal',num2str(n),'.jpeg'));

% imwrite(M1,filename,'jpeg')

axes(ax2);

axis tight

set(gca,'nextplot','replacechildren');

surf(sview');

shading interp;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j=13,k), time step = ',timestep],'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD');

xlabel('i coordinate');

ylabel('k coordinate');

set(gca,'fontname','Times New Roman','fontsize',10);

% F2 = getframe(gca,rect);

% M2 = frame2im(F2);

% filename = fullfile('C:\MATLAB7\work\cavity\figures',strcat('XZ_ETotal',num2str(n),'.jpeg')); % imwrite(M2,filename,'jpeg')

%***********************************************************************

% Cavity Power - Analitic expression

%***********************************************************************

axes(ax3);

% axis tight

% set(gca,'nextplot','replacechildren');

t = [1:1:param.nmax];

Psim = 1e3*Psim./max(Psim);

Panl = 1e3*Panl./max(Panl);

semilogy(t,Psim,'b.-',t,Panl,'rx-');

str(1) = {'Normalized Analytic Vs. Simulated Power'};

str(2) = {'as function of time-steps'};

title(str,'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD' );

xlabel('Time Step');

ylabel('Log(Power)');

legend('Simulation','Analytic','location','SouthEast');

set(gca,'fontname','Times New Roman','fontsize',10);

% F3 = getframe(gca,rectp);

% M3 = frame2im(F3);

% filename = fullfile('C:\MATLAB7\work\cavity\figures',strcat('Power',num2str(n),'.jpeg'));

% imwrite(M3,filename,'jpeg')

%文件3

% Source waveform function

function [source]=waveform(param,handles,n)

ax1 = handles.axes1;

type = get(handles.sourcetype,'value');

amp = get(handles.sourceamp,'value')*1e4;

f = get(handles.Frequency,'value')*1e9;

switch type

case 2

source = param.srcconst*(n-param.ndelay)*exp(-((n-param.ndelay)^2/param.tau^2));

case 1

source = param.srcconst*sin(param.dt*n*2*pi*f);

case 3

source = param.srcconst*exp(-((n-param.ndelay)^2/param.tau^2))*sin(2*pi*f*(n-param.ndelay)*param.dt) ;

otherwise

source = param.srcconst*(n-param.ndelay)*exp(-((n-param.ndelay)^2/param.tau^2)); end

%文件4

function [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,x,y,p)

a = param.a(x);

b = param.b(x);

c = param.c(y);

d = param.d(y);

hx(a+1:b,c:d,1:ke)=...

hx(a+1:b,c:d,1:ke)+...

param.db*(ey(a+1:b,c:d,2:kb)-...

ey(a+1:b,c:d,1:ke)+...

ez(a+1:b,c:d,1:ke)-...

ez(a+1:b,c+1:d+1,1:ke));

hy(a:b,c+1:d,1:ke)=...

hy(a:b,c+1:d,1:ke)+...

param.db*(ex(a:b,c+1:d,1:ke)-...

ex(a:b,c+1:d,2:kb)+...

ez(a+1:b+1,c+1:d,1:ke)-...

ez(a:b,c+1:d,1:ke));

hz(a:b,c:d,2:ke)=...

hz(a:b,c:d,2:ke)+...

param.db*(ex(a:b,c+1:d+1,2:ke)-...

ex(a:b,c:d,2:ke)+...

ey(a:b,c:d,2:ke)-...

ey(a+1:b+1,c:d,2:ke));

%文件5

function [ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,x,y,p)

a = param.a(x);

b = param.b(x);

c = param.c(y);

d = param.d(y);

if (param.is>=a)&(param.is<=b)|(param.js>=c)&(param.js<=d) source = waveform(param,handles,n);

else

source = 0;

end

ex(a:b,c+1:d,2:ke)=...

param.ca*ex(a:b,c+1:d,2:ke)+...

param.cb*(hz(a:b,c+1:d,2:ke)-...

hz(a:b,c:d-1,2:ke)+...

hy(a:b,c+1:d,1:ke-1)-...

hy(a:b,c+1:d,2:ke));

ey(a+1:b,c:d,2:ke)=...

param.ca*ey(a+1:b,c:d,2:ke)+...

param.cb*(hx(a+1:b,c:d,2:ke)-...

hx(a+1:b,c:d,1:ke-1)+...

hz(a:b-1,c:d,2:ke)-...

hz(a+1:b,c:d,2:ke));

ez(a+1:b,c+1:d,1:ke)=...

param.ca*ez(a+1:b,c+1:d,1:ke)+...

param.cb*(hx(a+1:b,c:d-1,1:ke)-...

hx(a+1:b,c+1:d,1:ke)+...

hy(a+1:b,c+1:d,1:ke)-...

hy(a:b-1,c+1:d,1:ke));

ez(param.is,param.js,1:ke)=ez(param.is,param.js,1:ke)+source; function FDTDonedimensionpipei(L,d,T)

%version1.0 终端匹配

%FDTDonedimensionpipei(6,0.18,0.5e-9)

t0=3*T;

c=3e8;

u=4*pi*1e-7;

e=8.8541878e-12;

dz=T*c/10;

Nz=L/dz;

Nz=fix(Nz);

dt=dz/2/c;

Ex=zeros(1,Nz+1);

B=zeros(1,Nz+1);

Hy=zeros(1,Nz);

Nt=2*Nz;

for n=0:Nt

t=n*dt;

F=exp(-(t-t0).^2./T^2);

Ex(1)=F;

for k=1:Nz

Hy(k)=Hy(k)+dt./u.*(Ex(k)-Ex(k+1))./dz;

end

for k=1:Nz-1

Ex(k+1)=Ex(k+1)+dt./e.*(Hy(k)-Hy(k+1))./dz;

end

Ex(1)=B(2)+(c*dt-dz)./(c*dt+dz).*(Ex(2)-B(1));

Ex(Nz+1)=B(Nz)+(c*dt-dz)./(c*dt+dz).*(Ex(Nz)-B(Nz+1));

Vref1=d.*Ex(Nz-300);

Vref2=d.*Ex(Nz-100);

plot(t,Vref1,'s');

hold on;

plot(t,Vref2,'rx');

hold on;

B=Ex;

end

function FDTD_debug

c_0 = 3E8; % Speed of light in free space

Nz = 100; % Number of cells in z-direction

Nt = 200; % Number of time steps

N = Nz; % Global number of cells

E = zeros(Nz,1); % E component

E2 = zeros(Nz,1);

E1 = zeros(Nz,1);

Es = zeros(Nz,Nz); % Es is the output for surf function

%H = zeros(Nz,1); % H component

pulse = 0; % Pulse

%to be set

mu_0 = 4.0*pi*1.0e-7; % Permeability of free space

eps_0 = 1.0/(c_0*c_0*mu_0); % Permittivty of free space

c_ref = c_0; % Reference velocity

eps_ref = eps_0;

mu_ref = mu_0;

f_0 = 10e9; % Excitation frequency

f_ref = f_0; % Reference frequency

omega_0 = 2.0*pi*f_0; % Excitation circular frequency

omega_ref = omega_0;

lambda_ref = c_ref/f_ref; % Reference wavelength

Dx_ref = lambda_ref/20; % Reference cells width

Dz = Dx_ref;

nDz = Dz/Dx_ref;

Z = Nz*Dz;

r = 1; % Normalized time step

Dt_ref = r*Dx_ref/c_ref; % Reference time step

Dt = Dt_ref;

% Source position

Nz_Source = 0.5*Nz;

N_Source = Nz_Source;

for n = 1:Nt-1

t = Dt_ref*r*(n-1); % Actual time

%Source function series

Source_type = 1;

switch Source_type

case 1 % modified source function

ncycles = 2;

if t < ncycles*2.0*pi/(omega_0)

pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t);

else

pulse = 0;

case 2 % sigle cos source function

if t < 2.0*pi/(omega_0)

pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);

else

pulse = 0;

end

case 3 % Gaussian pulse

if t < Dt_ref*r*50

pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2);

else

pulse = 0;

end

end

E2(N_Source) = E2(N_Source) - r*pulse;

E1 = E2;

E2 = E;

m = 2:Nz-1;

E(m) = r^2*(E2(m+1) - 2*E2(m) + E2(m-1)) + 2*E2(m) - E1(m);

%Boundary

E(1) = 0; E(Nz) = 0; % Dirichlet

%E(1) = E(2);E(Nz) = E(Nz-1); % Neumann

%E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz)); % Mur ABC right side z = Nz

%E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1)); % Mur ABC left side z = 0

%display

%choice***********************************************

display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf

switch display

case 1

i = 1:Nz;

plot(i, E(i));

axis([0 Nz -4 4]);

case 2

A = rot90(E);

imagesc(A);

case 3

i = 1:Nz;

for j = 1:Nz

Es(i,j) = E(i);

end

Es = rot90(Es);

j = 1:Nz;

surf(i,j,Es);

axis([0 Nz 0 Nz -4 4])

otherwise

A = rot90(E);

imagesc(A);

end

pause(0.005);

end

%***********************************************************************

% 3-D FDTD code with PEC boundaries

%*********************************************************************** %

% Program author: Susan C. Hagness

% Department of Electrical and Computer Engineering

% University of Wisconsin-Madison

% 1415 Engineering Drive

% Madison, WI 53706-1691

% 608-265-5739

% https://www.wendangku.net/doc/e39209758.html,

%

% Date of this version: February 2000

%

% This MATLAB M-file implements the finite-difference time-domain

% solution of Maxwell's curl equations over a three-dimensional

% Cartesian space lattice comprised of uniform cubic grid cells.

%

% To illustrate the algorithm, an air-filled rectangular cavity

% resonator is modeled. The length, width, and height of the

% cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and

% 2.0 cm (z-direction), respectively.

%

% The computational domain is truncated using PEC boundary

% conditions:

% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes

% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes

% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes

% These PEC boundaries form the outer lossless walls of the cavity.

%

% The cavity is excited by an additive current source oriented

% along the z-direction. The source waveform is a differentiated

% Gaussian pulse given by

% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),

% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-

% content pulse is approximately 7 GHz. The grid resolution

% (dx = 2 mm) was chosen to provide at least 10 samples per

% wavelength up through 15 GHz.

%

% To execute this M-file, type "fdtd3D" at the MATLAB prompt.

% This M-file displays the FDTD-computed Ez fields at every other

% time step, and records those frames in a movie matrix, M, which

% is played at the end of the simulation using the "movie" command.

%

%***********************************************************************

clear

%*********************************************************************** % Fundamental constants

%***********************************************************************

cc=2.99792458e8; %speed of light in free space

muz=4.0*pi*1.0e-7; %permeability of free space

epsz=1.0/(cc*cc*muz); %permittivity of free space

%*********************************************************************** % Grid parameters

%***********************************************************************

ie=50; %number of grid cells in x-direction

je=24; %number of grid cells in y-direction

ke=10; %number of grid cells in z-direction

ib=ie+1;

jb=je+1;

kb=ke+1;

is=26; %location of z-directed current source

js=13; %location of z-directed current source

kobs=5;

dx=0.002; %space increment of cubic lattice

dt=dx/(2.0*cc); %time step

nmax=500; %total number of time steps

%*********************************************************************** % Differentiated Gaussian pulse excitation

%***********************************************************************

rtau=50.0e-12;

tau=rtau/dt;

ndelay=3*tau;

srcconst=-dt*3.0e+11;

%*********************************************************************** % Material parameters

%***********************************************************************

eps=1.0;

sig=0.0;

%*********************************************************************** % Updating coefficients

%***********************************************************************

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));

cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));

da=1.0;

db=dt/muz/dx;

%*********************************************************************** % Field arrays

%***********************************************************************

ex=zeros(ie,jb,kb);

ey=zeros(ib,je,kb);

ez=zeros(ib,jb,ke);

hx=zeros(ib,je,ke);

hy=zeros(ie,jb,ke);

hz=zeros(ie,je,kb);

%*********************************************************************** % Movie initialization

%***********************************************************************

tview(:,:)=ez(:,:,kobs);

sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');

shading flat;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j,k=5), time step = 0']);

xlabel('i coordinate');

ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');

shading flat;

caxis([-1.0 1.0]);

MATLAB编程作业

《Matlab 编程训练》 作业 专 业 学生姓名 班级 学 号 指导教师 完成日期

实训一 MATLAB 语言介绍和数值计算 1．先求下列表达式的值，然后显示MATLAB 工作空间的使用情况并保存变量。 12 2sin851z e =+ . 2. 已知 1234413134787,2033657327A B --???? ????==???? ????-???? ，求下列表达式的值： (1) A+6*B 和A-B+I （其中I 为单位矩阵） A+6*B:

A-B+I: (2)A*B和A.*B A*B程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] c=A*B 结果： A.*B程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] D=A.*B 结果：

(3)A^3和A.^3 A^3程序： A=[12 34 -4;34 7 87;3 65 7] E=A^3 结果： A．^3程序： A=[12 34 -4;34 7 87;3 65 7] C=A.^3 (4)A/B及B\A A/B程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] C=A/B 结果：

B\A程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] D=B\A 结果： (5)将矩阵C=B\A的右下角2*2子矩阵赋给D, 并（3）保存变量（mat文件）程序： A=[12 34 -4;34 7 87;3 65 7]; B=[1 3 -1;2 0 3;3 -2 7]; C=B*inv(A); D=C(2:3,2:3) 结果：

MATLAB程序设计作业

Matlab程序设计 班级 姓名 学号

《MATLAB程序设计》作业 1、考虑如下x-y 一组实验数据： x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10] y=[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2] 分别绘出plot的原始数据、一次拟合曲线和三次拟合曲线，给出MATLAB代码和运行结果。 代码如下： x=[1,2,3,4,5,6,7,8,9,10]; y=[1.2,3,4,4,5,4.7,5,5.2,6,7.2]; plot(x,y); title('原始数据'); p=polyfit(x,y,1); q=polyval(p,x); figure,plot(x,q); title('一次拟合'); p=polyfit(x,y,2); q=polyval(p,x); figure,plot(x,q); title('二次拟合'); 运行结果如下：

1 2 3 4 5 6 7 8 9 10 12 3 4 5 6 7 8 原始数据 123 456789 102 2.5 3 3.54 4.5 55.56 6.57一次拟合 123456789 101 2 3 4 5 6 7 二次拟合 2、在[0，3π]区间，绘制y=sin(x)曲线（要求消去负半波，即(π，2π)区间内的函数值置零），求出曲线y 的平均值，以及y 的最大值及其最大值的位置。给出执行代码和运行结果。 代码如下： clear clc x=(0:0.01:3*pi); y=sin(x); plot(x,y); y1=(y>=0).*y; figure,plot(x,y1);

灰色预测模型的Matlab程序及检验程序(精)

灰色预测模型的Matlab 程序及检验程序 %灰色预测模型程序 clear syms a b; c=[a b]'; A=[46.2 32.6 26.7 23.0 20.0 18.9 17.5 16.3];% 原始序列 B=cumsum(A);%累加n=length(A); for i=1:(n-1) C(i)=(B(i)+B(i+1))/2; end %计算待定参数 D=A; D(1)=[]; D=D'; E=[-C; ones(1,n-1)]; c=inv(E*E')*E*D; c=c'; a=c(1); b=c(2); %预测往后预测5个数据 F=[];F(1)=A(1); for i=2:(n+5) F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a; end G=[];G(1)=A(1); for i=2:(n+5) G(i)=F(i)-F(i-1); end t1=2002:2009; t2=2002:2014; G plot(t1,A,'o',t2,G) %灰色预测模型检验程序 function [ q,c,p ] = checkgm( x0,x1 ) %GM 检验函数 %x0 原始序列

%x1 预测序列 %·返回值 % q –- 相对误差 % c -- ·方差比 % p -- 小误差概率 e0=x0-x1; q=e0/x0; s1=var(x0); %qpa=mean(e0); s2=var(e0); c=s2/s1; len=length(e0); p=0; for i=1:len if(abs(e0(i)) < 0.6745*s1) p=p+1; end end p=p/len; end

MATLAB第一章作业答案

第一章 M A T L A B 概况与基本操作 1.选择题（每题2分，共20分）： (1)最初的MATLAB 核心程序是采用D 语言编写的。 A.PASCAL B.C C.BASIC D.FORTRAN (2)即将于2011年9月发布的MATLAB 新版本的编号为C 。 A.MATLAB 2011Ra B.MATLAB 2011Rb C.MATLAB R2011b D.MATLAB R2011a (3)在默认设置中，MATLAB 中的注释语句显示的颜色是B 。 A.黑色 B.绿色 C.红色 D.蓝色 (4)如果要以科学计数法显示15位有效数字，使用的命令是A 。 A.format long e B.format long C.format long g D.format long d (5)在命令窗口新建变量a 、b ，如果只查看变量a 的详细信息，使用的命令为A 。 A.whos a B.who a C.who D.whos (6)如果要清除工作空间的所有变量，使用的命令为C 。 A.clear B.clear all C.两者都可 D.两者都不可 (7)在创建变量时，如果不想立即在命令窗口中输出结果，可以在命令后加上B 。 A.冒号 B.分号 C.空格 D.逗号 (8)如果要重新执行以前输入的命令，可以使用D 键。 A.下箭头↓ B.右箭头→ C.左箭头← D.上箭头↑ (9)如果要查询函数det 的功能和用法，并显示在命令窗口，应使用命令C 。 A.doc B.lookfor C.help D.三者均可 (10)如果要启动Notebook 文档，下列D 操作是可行的。 A.在命令窗口输入notebook 命令 B.在命令窗口输入notebook filename 命令 C.在Word 中启动M-book 文档 D.三者均可 2.填空题（每空1分，共20分）： (1)MATLAB 是matrix 和laboratory 两个单词前三个字母的组合，意为“矩阵实验室”，它的创始人是Cleve Moler 和Jack Little 。 (2)在MATLAB 的默认设置中，关键字显示的字体为蓝色，命令、表达式、计算结果显示的字体为黑色，字符串显示的字体为褐红色，注释显示的字体为绿色，错误信息显示的字体为红色。 (3)在命令窗口中，输出结果显示为各行之间添加空行的命令为format loose ，各行之间不添加空行的命令为format compact 。 (4)在MATLAB 中，各种标点符号的作用是不同的。例如，空格的作用是分隔数组每行各个元素，逗号的作用是分隔数组每行各个元素或函数的各个输入参数，分号的作用是作为不显示命令结果的命令行的结尾或分隔数组各列，冒号的作用是生成一维数组或表示数组全部元素，百分号的作用是引导一行注释，…的作用是连接相邻两行，感叹号的作用是调用操作系统命令。 3.程序设计题（每题10分，共40分） (1)以25m/s 的初速度向正上方投球（g=9.8m/s 2），计算到达最高点的时间tp 以及球从出发点到 最高点的距离hp 。 解：根据物理学知识，物体上抛运动的速度与经过的时间之间的关系为0p p v v gt =-，因此所需要的时间为0p p v v t g -=。而到达最高点时的速度0p v =，因此可根据此公式求出tp ： v0=25;g=9.8;vp=0; tp=(v0-vp)/g tp = 2.5510

灰色预测MATLAB程序

作用：求累加数列、求a b的值、求预测方程、求残差 clc %清屏，以使结果独立显示 x=[ ]; format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x'; end n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值，并将值赋予矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end

for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值 for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值 ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,: );%显示输出预测值的累加数列 end var(1,:)=ago(1,:) %显示输出预测值 for i=1:n %如改n为n+m-1，可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值end for i=1:n error(i,:)=x(i,:)-var(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列 alpha %显示输出参数数列 var %显示输出预测值 error %显示输出误差 c %显示后验差的比值 作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程 clc,clear x0=[ ]'; %注意这里为列向量 n=length(x0); lamda=x0(1:n-1)./x0(2:n) %计算级比 range=minmax(lamda') %计算级比的范围 x1=cumsum(x0) %累加运算 B=[*(x1(1:n-1)+x1(2:n)),ones(n-1,1)]; Y=x0(2:n); u=B\Y %拟合参数u(1)=a,u(2)=b x=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解

灰色预测matlab程序

Matlab的灰色预测程序： y=input('请输入数据'); n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i) end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1))/2; B(i,2)=1; end BT=B'; for j=1:(n-1) YN(j)=y(j+1); end YN=YN'; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; t_test=input('输入需要预测的个数'); i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test; yn=ys(2:n+t_test); plot(x,y,'^r',xs,yn,'*-b'); det=0; for i=2:n det=det+abs(yn(i)-y(i)); end det=det/(n-1); disp(['百分绝对误差为：',num2str(det),'%']); disp(['预测值为：',num2str(ys(n+1:n+t_test))]); 请输入数据[29.8 30.11 41.05 70.12 77.79 77.79 104.82 65.22 82.7 100.79] 输入需要预测的个数4

百分绝对误差为：14.5128% 预测值为：110.5718 120.8171 132.0116 144.2434

自动控制原理Matlab程序作业(精)

自控控制原理 MATLAB 程序设计作业 指导老师:汪晓宁 目录 一、题目 (2) 二、运行结果 (3) 三、程序说明 (8) 四、附录 ............................................ 9 代码 . ............................................. 9 参考文献 .. (17) 一、题目 用 Matlab 创建用户界面,并完成以下功能 a 将产生未综合系统的根轨迹图以及 0.707阻尼比线, 你可以交互地选择交点的运行点。界面能显示运行点的坐标、增益值以及近似为二阶系统估算的超调量、调整时间、峰值时间、阻尼比、无阻尼自然震荡频率以及稳态误差 b 显示未综合系统的阶跃响应 c 输入控制器的参数, 绘制综合后系统的根轨迹图以及显示综合的设计点 (主导极点 , 允许不断改变控制器参数,知道所绘制的根轨迹通过设计点 d 对于综合后的系统, 显示运行点的坐标、增益,近似为二阶系统估算的超调量、调整时间、峰值时间、阻尼比、无阻尼自然震荡频率以及误差系数 e 显示综合后系统的阶跃响应 二、运行结果

输入传递函数分子分母 生成根轨迹图

选择点并得到该点各项参数在下方输出面板输出 获得阶跃响应图 用 rltool(辅助,选择合适的插入零点

输入零点,并得到根轨迹图

选择根轨迹图上的任一点,得到数据,在下方输出面板输出得到阶跃响应图 三、运行说明

第一步, 在请输入分子后的输入框输入传递函数分子的矩阵, 在下一输入框输入传递函数分母并按“生成根轨迹图”按钮获得根轨迹 第二步, 按选择点并显示各参数获得根轨迹图上任一点的各项数据, 数据全部输出在下方输出面板 第三步,按“生成阶跃响应图”按钮可以获得该函数的阶跃响应 第四步,在“请输入插入零点”后的输入框中输入参数,并按“生成综合后根轨迹图” 按钮产生根轨迹 (可以通过点击“根轨迹校正”按钮,调用工具箱拖动零点进行快速查看根轨迹图,选择合适的根轨迹再在输入框中输入零点的值 第五步,按“选择点并显示各参数(综合后系统”选取各点,查阅参数,数据输出在下方输出面板上 第六步,按“生成阶跃响应图(综合后系统”可以得到综合后系统的阶跃响应 最后,点击“退出”结束程序 四、附录 代码: function varargout = Liushuai20122510(varargin % LIUSHUAI20122510 MATLAB code for Liushuai20122510.fig % LIUSHUAI20122510, by itself, creates a new LIUSHUAI20122510 or raises the existing % singleton*. %

灰色预测模型matlab程序精确版

灰色预测模型matlab程序 %下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,mat lab6.5 ,使用本程序请注明，程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN; au = au0'; %B,au0,au

ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '('; rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,r ightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1; au2 = au20'; %z4,X1,G,au20

Matlab程序设计大作业(终审稿)

M a t l a b程序设计大作 业 公司内部档案编码：[OPPTR-OPPT28-OPPTL98-OPPNN08]

Matlab程序设计 课程大作业 题目名称：_________________________________ 班级：_________________________________ 姓名：_________________________________ 学号：_________________________________ 课程教师：温海骏 学期： 2015-2016学年第2学期 完成时间：

MATLAB优化应用 §1 线性规划模型 一、线性规划问题： 问题1：生产计划问题 假设某厂计划生产甲、乙两种产品，现库存主要材料有A类3600公斤，B类2000公斤，C类3000公斤。每件甲产品需用材料A类9公斤，B类4公斤，C类3公斤。每件乙产品，需用材料A类4公斤，B类5公斤，C类10公斤。甲单位产品的利润70元，乙单位产品的利润120元。问如何安排生产，才能使该厂所获的利润最大。 问题2：投资问题 某公司有一批资金用于4个工程项目的投资，其投资各项目时所得的净收益(投入资金百分比)如下表：工程项目收益表 由于某种原因，决定用于项目A的投资不大于其他各项投资之和而用于项目B和C的投资要大于项目D的投资。试确定该公司收益最大的投资分配方案。 问题3：运输问题

有A 、B 、C 三个食品加工厂，负责供给甲、乙、丙、丁四个市场。三个厂每天生产食品箱数上限如下表： 四个市场每天的需求量如下表： 从各厂运到各市场的运输费(元/每箱)由下表给出： 求在基本满足供需平衡的约束条件下使总运输费用最小。 §2 多目标规划模型 多目标规划定义为在一组约束下，多个不同的目标函数进行优化设计。 数学模型： 12min ()() ().()0,1,2, ,m j f x f x f x st g x j k ???? ≤= 其中x=(x 1 ,x 2 , … ,x n )为一个n 维向量；f i (x)为目标函数，i=1, 2, … ,m; g j (x)为系统约束, j=1, 2, … ,k 。

灰色预测模型matlab程序精确版

%x=[1019,1088,1324,1408,1601];gm1(x); 测试数据%二次拟合预测GM(1,1)模型 function gmcal=gm1(x) if nargin==0 x=[1019,1088,1324,1408,1601] end format long g sizex=length(x); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=length(z1); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN; au = au0'; %B,au0,au afor = au(1); ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '(';

rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+ ',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1=length(ze1); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1; au2 = au20'; %z4,X1,G,au20 Aval = au2(1); Bval = au2(2); %Aval,Bval %输出预测的 A,B的值 strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',lef tbra,num2str(Bval),rightbra) %输出时间响应方程 nfinal = sizex-1 + 1;（其中+1可改为+5等其他数字，即可预测更多的数字） %决定预测的步骤数5 这个步骤可以通过函数传入 %nfinal = sizexd2 - 1 + 1; %预测的步骤数 1 for k3=1:nfinal x3fcast(k3) = constant1*exp(afor1*k3)+ua; end %x3fcast %一次拟合累加值 for k31=nfinal:-1:0 if k31>1 x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); else if k31>0

matlab程序设计作业

Matlab程序设计作业 姓名: 学号: 专业:

? MATLAB 程序设计》作业 1、考虑如下x-y 一组实验数据: x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10] y 二[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2] 分别绘出plot 的原始数据、一次拟合曲线和三次拟合曲线，给出 原始曲线 MATLAB 代码和运行结果。 7 6 5 4 3 2 2 3 4 5 6 7 8 9 10

7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 10 一次拟合 三次拟合

x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; y=[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2]; figure; plot(x,y) p1=polyfit(x,y,1); y1=polyval(p1,x); figure; plot(x,y1) p2=polyfit(x,y,3); y2=polyval(p2,x); figure; plot(x,y2) 2、在［0, 3n区间，绘制y二Sin(x)曲线(要求消去负半波，即(n 2n)区间内的函数值置零)，求出曲线y 的平均值，以及y 的最大值及其最大值的位置。给出执行代码和运行结果。 x=0:pi/1000:3*pi; y=Sin(x); y1=(y>=0).*y; %消去负半波figure(1); plot(x,y1, 'b' ); a=mean(y1) %求出y1 的平均值 b=max(y1) %求出y1 的最大值b, 以及最大值在矩阵中的位置; d=x(find(y1==b)) >> ex1 a = 0.4243 b = 1 d = 1.5708 7.8540 >>

MATLAB第一章作业答案

第一章M A T L A B概况与基本操作 1.选择题： (1)最初的MATLAB核心程序是采用A语言编写的。 A.FORTRAN B.C C.BASIC D.PASCAL (2)即将于2011年9月发布的MATLAB新版本的编号为D。 A.MATLAB 2011Ra B.MATLAB 2011Rb C.MATLAB R2011a D.MATLAB R2011b (3)在默认设置中，MATLAB中的注释语句显示的颜色是D。 A.黑色 B.蓝色 C.红色 D.绿色 (4)如果要以科学计数法显示15位有效数字，使用的命令是B。 A.format long B.format long e C.format long g D.format long d (5)在命令窗口新建变量a、b，如果只查看变量a的详细信息，使用的命令为B。 A.who a B.whos a C.who D.whos (6)如果要清除工作空间的所有变量，使用的命令为 C 。 A.clear B.clear all C.两者都可 D.两者都不可 (7)在创建变量时，如果不想立即在命令窗口中输出结果，可以在命令后加上D。 A.冒号 B.逗号 C.空格 D.分号 (8)如果要重新执行以前输入的命令，可以使用B键。 A.下箭头↓ B.上箭头↑ C.左箭头← D.右箭头→ (9)如果要查询函数inv的相关信息，并显示在命令窗口，应使用命令A。 A.help B.lookfor C.doc D.三者均可 (10)如果要启动Notebook文档，下列D操作是可行的。 A.在命令窗口输入notebook命令 B.在命令窗口输入notebook filename命令 C.在Word中启动M-book文档 D.三者均可 2.填空题： (1)MATLAB是MATrix和LABoratory两个单词前三个字母的组合，意为“矩阵实验室”，它的创始人是Cleve Moler和Jack Little。 (2)在MATLAB的默认设置中，关键字显示的字体为蓝色，命令、表达式、计算结果显示的字体为黑色，字符串显示的字体为紫色，注释显示的字体为绿色，错误信息显示的字体为红色。 (3)在命令窗口中，输出结果显示为各行之间添加空行的命令为format compact，各行之间不添加空行的命令为format compact。 备注：本题布置给大家时有一点小错误，现在予以更正。 (4)在MATLAB中，各种标点符号的作用是不同的。例如，空格的作用是分隔数组中每一行的各个元素，逗号的作用是分隔数组中每一行的各个元素或不同的命令，分号的作用是分隔数组中的各行或控制命令执行结果是否在命令窗口显示，冒号的作用是生成一维数组或显示全部元素，百分号的作用是注释行的开头，…的作用是把相邻两行的语句连接为一个命令，感叹号的作用是执行操作系统命令。 3.先建立自己的工作目录，再将自己的工作目录设置到MATLAB搜索路径下。请写出操作步骤或用Matlab命令实现。用help命令能查询到自己的工作目录吗？ 解：操作步骤： (1)在Windows环境中建立一个工作目录，如：c:\mywork； (2)启动MATLAB，在File菜单中选择Set Path…命令，显示出如下图的对话框：

灰色预测MATLAB程序

灰色预测 专设工⑼他QA—叫吋)为原始数列.其1次累 ?加生成数列为恥=妙①曲⑵，…卅何),其中 X° 仇)二工* ° (0.址=1=2= -:n 5-1 卷定义卫的决导数为 d(k) = *町(上)=x 叫咼-x cl)(Jt-l). 令为数列工①的邻值生成数列.即 却(去)=^(*) + (1- a)x0)(t-lX 于是定义GM (L 1)的灰微分方程模型为 d(k)-血⑴住)=K 即或严＞(￡) + “尹⑻=人⑴ 在式(1)中』。＞(灼称为灰导数，我称为发展系数, 弧称为白化背景值，b称为灰作用量乜 将时刻表殳二2「3「/代入(1)式有 V!1「— a y=代⑶ B =I b*- : X闵0)-Z，：］(K)1

于是G\I <1?1）複至可表示为Y = Bu. 現在问题归结为求sb 在值。用一元线性回归?即最小二秦法求它们的活计值 为 注二实陌上回归分析中求估计值是用软件计尊的?有标准程序求解，iOmaClab 等。 GM <1? 1>的白化晏 対于G\I <1> 1）的灰微分方程（1） >如果将灰导数打（Q 的时刻 视为连绫变里"则x°）视为时问（函数卅⑺,于是*〉（Q 対血于导数里级 心2 >白化背臬值申的对应于导数卅⑴。于是G\I （1，1）的坝徽 分方樂対应于的白微分 方程为 内?则数堀列X?可以塗互G\I <19 1） 且可以进行页 色预测。否朋，対数摄做适当的克换处理■如平移叢换: 取C 使得鞍据列 严伙）=工⑴伙）+ G 上=1,2,…, 的级比都華住可吝禎盖内。 心⑴⑴ + o?i> (r)二dr <2) GM mi ）质色预测的步骤 1 ?教摇的枪绘与处連 为了ftilGAl （1,1）建複方法的可行性，亲要为已知期S 做必要的检蛉 处理。 设原始教据列为了 逛=（乂°（1）*6（2）严炉00； >计算数列的级比 如果所有的级比都落在可容覆盖区间 ? fc = A- 2,3"?

Matlab作业习题与答案详解(附程序)

clear all;clc;close all; x=-10:0.01:20; y=4*sin(x)./x; ymin=min(y) 二、蒙特卡罗算法的数值计算 当前的油位高度是2.3米，见图1。模拟油流进储油罐的过程（图维数任选），请计算罐内油量。三维的效果图参见图2。储油罐由两部分组成，中间是圆柱体，两端是球罐体。（本题简化自2011年UCMCM A题《储油罐的变位识别与罐容表标定》，细节参见原题原题附件2 cumcm2010A.doc。） 图1

图2 主程序： clc; clear all; close all; center1=[-3.375,0,1.5]; %左球罐中心center2=[3.375,0,1.5]; %右球罐中心 n=10000; %每次的撒点数 delta=0.02; %层高 h=3; en=h/delta; Show; %画出油罐 for i=0:en-1 x=(rand(1,n)-0.5)*10; %随机生成n个点

y=(rand(1,n)-0.5)*h; z=(rand(1,n)*delta+i*delta); Z=[x;y;z]; [dis1 dis2]=juli(center1,center2,Z); %算出各点对应的距离 index=find(((x>-4&x<4)&dis2<1.5)|(x<-4|x>4)&dis1<1.625); %找出在罐内的点 plot3(x(index),y(index),z(index),'.k'); %画出在罐内的点 drawnow end 子程序1： function [dis1 dis2]=juli(a,b,q) d11=q(1,:)-a(1); d12=q(2,:)-a(2); d13=q(3,:)-a(3); d1=sqrt(d11.^2+d12.^2+d13.^2); d21=q(1,:)-b(1); d22=q(2,:)-b(2);

灰色预测MATLAB程序

灰色预测M A T L A B程 序 文档编制序号：[KKIDT-LLE0828-LLETD298-POI08]

灰色预测 作用：求累加数列、求a b的值、求预测方程、求残差 clc %清屏，以使结果独立显示 x=[ ]; format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x'; end n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值，并将值赋予矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值 for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值 ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i- 1))+alpha(2,:)/alpha(1,:);%显示输出预测值的累加数列 end var(1,:)=ago(1,:) %显示输出预测值 for i=1:n %如改n为n+m-1，可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值end for i=1:n error(i,:)=x(i,:)-var(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列

Matlab振动程序-代码作业

一、课题任务要求 随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。 本文为机械测试信号的时域和频域分析，其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。在进行上述分析之前先要对振动信号进行拟合。机械振动分为确定性振动和随机振动，确定性振动又分为周期振动和非周期振动，周期振动又进一步分为简谐振动和复杂的周期振动。所以可以根据上述的分类来拟合振动信号。在设计信号的处理程序时，用MATLAB中的相关函数来对所拟合的振动信号进行时域分析和频域分析，并对绘出的频谱图进行说明。 二、技术路线 对机械振动信号的时域和频域采集，根据所拟合的振动信号，选取所需要的时域性能指标和频域分析的性能指标对振动信号进行分析。其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。 现构造一个振动信号（在该程序中以两个衰减振动分量和一个随

机数rand 之和来拟合振动信号），再利用MATLAB 中的函数mean （）、min （）、max （）、std （）对离散序列中的平均值、最大值、最小值、标准差等时域性能进行分析，通过调用函数fft （y ）；psd （y ）；rcep （y ）对该振动信号进行频域内的性能分析。 在设计过程中的理论知识有离散傅立叶变换（DFT ）、功率谱的概念和意义以及倒谱的概念和意义。 ①DFT 的定义和意义： DFT 的定义式为： DFT 的意义：DFT 的意义在于它表示信号中的各个频率的分量的振动幅值的大小，亦即该分量对于振动信号影响的大小。通过DFT 的快速算法FFT 可以很方便的将振动信号的各个分量的幅值比重计算出来。即可以把信号的主频分量提取出来。 ②功率谱的概念和意义： 功率谱的定义式为:若X (Ω)=DFT[x(m)]，x(n)为N 点序列。则 X * （Ω） =DFT[x N (-m)] =0k

灰色预测模型的MATLAB 程序及检验程序

灰色预测模型的Matlab程序及检验程序%灰色预测模型程序 clear syms a b; c=[a b]'; A=[46.232.626.723.020.018.917.516.3];%原始序列B=cumsum(A);%累加 n=length(A); for i=1:(n-1) C(i)=(B(i)+B(i+1))/2; end %计算待定参数 D=A; D(1)=[]; D=D'; E=[-C;ones(1,n-1)]; c=inv(E*E')*E*D; c=c'; a=c(1); b=c(2); %预测往后预测5个数据 F=[];F(1)=A(1); for i=2:(n+5) F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a; end G=[];G(1)=A(1); for i=2:(n+5) G(i)=F(i)-F(i-1); end t1=2002:2009; t2=2002:2014; G plot(t1,A,'o',t2,G) %灰色预测模型检验程序 function[q,c,p]=checkgm(x0,x1) %GM检验函数 %x0原始序列 %x1预测序列 %·返回值

%q–-相对误差 %c--·方差比 %p--小误差概率 e0=x0-x1; q=e0/x0; s1=var(x0); %qpa=mean(e0); s2=var(e0); c=s2/s1; len=length(e0); p=0; for i=1:len if(abs(e0(i))<0.6745*s1) p=p+1; end end p=p/len; end 等级相对误差q方差比C小误差概论P I级<0.01<0.35>0.95 II级<0.05<0.50<0.80 III级<0.10<0.65<0.70 IV级>0.20>0.80<0.60

电气工程软件训练(三)——Matlab 作业

D1( 江苏大学 《电气工程软件训练三》课程设计报告 设计题目：MATLAB 专业班级：J电气1401 学生姓名：唐鹏 学生学号：4141127007 指导老师： 完成日期： 江苏大学京江学院

一MATLAB课程设计的目的和要求 1.MATLAB软件功能简介 MATLAB的名称源自Matrix Laboratory，1984年由美国Mathworks公司推向市场。它是一种科学计算软件，专门以矩阵的形式处理数据。MATLAB将高性能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛的应用于科学计算、控制系统和信息处理等领域的分析、仿真和设计工作。MATLAB 软件包括五大通用功能，数值计算功能（Nemeric）、符号运算功能（Symbolic）、数据可视化功能（Graphic）、数字图形文字统一处理功能（Notebook）和建模仿真可视化功能（Simulink）。其中，符号运算功能的实现是通过请求MAPLE内核计算并将结果返回到MATLAB命令窗口。该软件有三大特点，一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks公司已推出30多个应用工具箱。MATLAB在线性代数、矩阵分析、数值及优化、数值统计和随机信号分析、电路与系统、系统动力学、次那好和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。 2.MATLAB课程设计的目的 本次课程设计主要是为了使学生了解MATLAB软件的基本知识，熟悉MATLAB的上机环境，掌握MATLAB数值运算、程序设计、二维/三维绘图、符号运算、Simulink仿真等相关知识，并初步具备将一般数学问题转化为对应的计算机进行处理的能力，以便为今后进一步的学习打下坚定基础。 二MATLAB课程内容 1 MATLAB语言基础 实验目的：基本掌握MATLAB 向量、矩阵、数组的生成及其基本运算(区分数组运算和矩阵运算)、常用的数学函数。了解字符串的操作。 实验内容: ①创建以下矩阵：A为初值为1，终值为12，元素数目为6的行向量；

灰色预测模型matlab程序精确版

灰色预测模型matlab程序 灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法，其预测数据与原始数据存在一定的规律相似性 %下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,mat lab6.5 ,使用本程序请注明，程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN;

%B,au0,au afor = au(1); ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '('; rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,r ightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1;