搜档网
当前位置:搜档网 › 数字信号处理实验全部程序MATLAB

数字信号处理实验全部程序MATLAB

数字信号处理实验全部程序MATLAB
数字信号处理实验全部程序MATLAB

实验一熟悉MATLAB环境

一、实验目的

(1)熟悉MATLAB的主要操作命令。

(2)学会简单的矩阵输入和数据读写。

(3)掌握简单的绘图命令。

(4)用MATLAB编程并学会创建函数。

(5)观察离散系统的频率响应。

二、实验内容

认真阅读本章附录,在MATLAB环境下重新做一遍附录中的例子,体会各条命令的含义。在熟悉了MATLAB基本命令的基础上,完成以下实验。

上机实验内容:

(1)数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4],B=[3 4 5 6],求C=A+B,D=A-B,E=A.*B,F=A./B,G=A.^B并用stem语句画出A、B、C、D、E、F、G。

实验程序:

A=[1 2 3 4];

B=[3 4 5 6];

n=1:4;

C=A+B;D=A-B;E=A.*B;F=A./B;G=A.^B;

subplot(4,2,1);stem(n,A,'fill');xlabel ('时间序列n');ylabel('A');

subplot(4,2,2);stem(n,B,'fill');xlabel ('时间序列n ');ylabel('B'); subplot(4,2,3);stem(n,C,'fill');xlabel ('时间序列n ');ylabel('A+B'); subplot(4,2,4);stem(n,D,'fill');xlabel ('时间序列n ');ylabel('A-B'); subplot(4,2,5);stem(n,E,'fill');xlabel ('时间序列n ');ylabel('A.*B'); subplot(4,2,6);stem(n,F,'fill');xlabel ('时间序列n ');ylabel('A./B'); subplot(4,2,7);stem(n,G,'fill');xlabel ('时间序列n ');ylabel('A.^B');

运行结果:

(2)用MATLAB实现以下序列。

a)x(n)=0.8n0≤n≤15

实验程序:

n=0:15;x=0.8.^n;

stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x(n)=0.8^n');

b)x(n)=e(0.2+3j)n0≤n≤15

实验程序:

n=0:15;x=exp((0.2+3*j)*n);

stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x(n)=exp((0.2+3*j)*n)'); 运行结果:

a)的时间序列 b)的时间序列c)x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π) 0≤n≤15

实验程序:

n=0:1:15;

x=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi);

stem(n,x,'fill'); xlabel('时间序列n ');

ylabel('x(n)=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi)');

运行结果:

d)将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期

实验程序:

n=0:1:63;

x=3*cos(0.125*pi*rem(n,16)+0.2*pi)+2*sin(0.25*pi*rem(n,16)+0.1*pi); stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x16(n)');

e)将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期

实验程序:

n=0:1:39;

x=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.25*pi*rem(n,10)+0.1*pi); stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x10(n)');

运行结果:

d )的时间序列

e )的时间序列

(3)x(n)=[1,-1,3,5],产生并绘出下列序列的样本。

a )x 1(n)=2x(n+2)-x(n-1)-2x(n)

实验程序:

n=0:3;

x=[1 -1 3 5];

x1=circshift(x,[0 -2]);x2=circshift(x,[0 1]);x3=2*x1-x2-2*x;

stem(x3,'fill'); xlabel ('时间序列n

');ylabel('x1(n)=2x(n+2)-x(n-1)-2x(n)');

b )∑=-=5

1k 2)k n (nx (n) x

实验程序:

n=0:3;

x=[1 -1 3 5];

x1=circshift(x,[0 1]);x2=circshift(x,[0 2]);x3=circshift(x,[0 3]);

x4=circshift(x,[0 4]);x5=circshift(x,[0 5]);

xn=1*x1+2*x2+3*x3+4*x4+5*x5;

stem(xn,'fill'); xlabel('时间序列n ');

ylabel('x2(n)=x(n-1)+2x(n-2)+3x(n-3)+4x(n-4)+5x(n-5)');

运行结果:

a)的时间序列b)的时间序列(4)绘出时间函数的图形,对x轴、y轴图形上方均须加上适当的标注。

a) x(t)=sin(2πt) 0≤t≤10s b) x(t)=cos(100πt)sin(πt) 0≤t≤4s

实验程序:

clc;

t1=0:0.001:10;t2=0:0.01:4;

xa=sin(2*pi*t1);xb=cos(100*pi*t2).*sin(pi*t2);

subplot(2,1,1);

plot(t1,xa);xlabel ('t');ylabel('x(t)');title('x(t)=sin(2*pi*t)'); subplot(2,1,2);

plot(t2,xb);xlabel

('t');ylabel('x(t)');title('x(t)=cos(100*pi*t2).*sin(pi*t2)');

运行结果:

(5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1

实验程序:

clc;

n1=input('请输入起点:');

n2=input('请输入终点');

n0=input('请输入阶跃位置');

n=n1:n2;

x=[n-n0>=0];

stem(n,x,'fill');xlable('时间序列n');ylable('u(n-n0)');

请输入起点:2

请输入终点:8

请输入阶跃位置:6

运行结果:

(5)运行结果 (6)运行结果

(6)给一定因果系统)0.9z 0.67z -1)/(1z 2(1H(z)-2-1-1+++

=求出并绘制H(z)的幅

频响应与相频响应。

实验程序:

a=[1 -0.67 0.9];

b=[1 sqrt(2) 1];

[h w]=freqz(b,a);

fp=20*log(abs(h));

subplot(2,1,1); plot(w,fp);xlabel('时间序列t');ylabel('幅频特性');

xp=angle(h);

subplot(2,1,2);

plot(w,xp);xlabel('时间序列t');ylabel('相频特性');

运行结果:(右上图)

常用典型序列

单位采样序列

function [x,n]=impseq(n1,n2,n0)

n=[n1:n2];

x=[(n-n0)==0];

[x,n]=impseq(-2,8,2);

stem(n,x);

title('电信1201')

n0=-2;

n=[-10:10];

nc=length(n);

x=zeros(1,nc);

for i=1:nc

if n(i)==n0

x(i)=1

end

end

stem(n,x) ;

title('电信1201 采样序列第二种方法')

单位阶跃序列

function [x,n]=stepseq(n1,n2,n0) n=[n1:n2];

x=[(n-n0)>=0];

[x,n]=stepseq(-2,8,2);

stem(n,x);

title('电信1201')

实数指数

n=[0:10];

x=0.9.^n;

stem(n,x);

title('电信1201')

复数指数序列

n=[-10:10];

alpha=-0.1+0.3*j;

x=exp(alpha*n);

real_x=real(x); image_x=imag(x);

mag_x=abs(x); phase_x=angle(x);

subplot(2,2,1); stem(n,real_x) ;title('电信1201')

subplot(2,2,2); stem(n,image_x);title('电信1201')

subplot(2,2,3); stem(n,mag_x);title('电信1201')

subplot(2,2,4); stem(n,phase_x);title('电信1201')

正余弦

n=[0:10];

x=3*cos(0.1*pi*n+pi/3);

stem(n,x);

title('电信1201')

例:求出下列波形

x1(n)=2x(n-5)-3x(n+4)

function [y,n]=sigadd(x1,n1,x2,n2)

n=min(min(n1),min(n2)) : max(max(n1),max(n2));

y1=zeros(1,length(n));

y2=y1;

y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;

y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;

y=y1+y2;

function [y,n]=sigshift(x,m,n0)

n=m+n0;

y=x;

n=[-2:10];

x=[1:7,6:-1:1];

[x11,n11]=sigshift(x,n,5);

[x12,n12]=sigshift(x,n,-4);

[x1,n1]=sigadd(2*x11,n11,-3*x12,n12);

stem(n1,x1);

title('电信1201')

已知某一系统方程为:y[n]-y[n-1]+0.9y[n-2]=x[n]计算并画出脉冲响应h(n),n=(-20,100) n=(-20:100);

num=(1); den=[1 -1 0.9];

x=impseq(-20,100,0);

h=filter(num,den,x);

stem(n,h)

xlabel('时间序号N');

ylabel('脉冲响应h');

title('脉冲响应电信1201');

实现两个序列a,b的卷积

x=[3,11,7,0,-1,4,2];

h=[2,3,0,-5,2,1];

c=conv(x,h);

stem(c);

title('电信1201')

自己给的数

x=[0.5,1,1.5,0,0];

h=[1,1,1];

c=conv(x,h);

stem(c);

title('电信1201')

试求卷积C(t)=f1(t)*f2(t),并绘制出f1、f2、及卷积以后的波形。function [y,ny]= conv_m(x,nx,h,nh,p) % 信号处理的改进卷积程序

nyb=nx(1)+nh(1);

nyc=nx(length(x))+nh(length(h));

ny=(nyb :p: nyc);

y=conv(x , h);

p=0.1;

t1= [0:p:1];

f1=t1.*(t1>0);

t2= [-1:0.1:2];

f2=t2.*exp(t2).*(t2>=0)+exp(t2).*(t2<0);

[y,ny]=conv_m(f1,t1,f2,t2,p);

subplot(3,1,1);stem(t1,f1);title('电信1201')

subplot(3,1,2);stem(t2,f2);title('电信1201')

subplot(3,1,3); stem(ny,y);title('电信1201')

实验二

function xk=dfs(xn,N)

n=(0:1:N-1);

k=n;

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^nk;

xk=xn* WNnk;

xn=[0,1,2,3];

N=4;

xk=dfs(xn,N)'

IDFS

function xn=idfs(xk,N)

n=(0:1:N-1);

k=n;

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^(-nk);

xn=xk*WNnk/N;

分析:因为x(n)是复指数,它满足周期性,我们将在两个周期中的401个频点上作计算来观n=0:10;

x=(0.9*exp(j*pi/3)).^n;

k=-200:200;

w=(pi/100)*k;

X=x*(exp(-j*pi/100)).^(n'*k);

magX=abs(X);

angX=angle(X);

subplot(2,1,1);

plot(w/pi,magX);

subplot(2,1,2);

plot(w/pi,angX/pi);

title('电信1201')

检验频移特性

n=0:100; x=cos(pi*n/2);

k=-100:100; w=(pi/100)*k;

X=x*(exp(-j*pi/100)).^(n'*k);

y=exp(j*pi*n/4).*x;

Y=y*(exp(-j*pi/100)).^(n'*k);

subplot(2,2,1);plot(w/pi,abs(X));

axis([-1,1,0,60]);

title('幅值电信1201');xlabel('以pi为单位的频率');ylabel('绝对值X');

subplot(2,2,2);plot(w/pi,angle(X)/pi);

title('相位电信1201');xlabel('以pi为单位的频率');ylabel('绝对值X');

axis([-1,1,-1,1]);

subplot(2,2,3);plot(w/pi,abs(Y));

title('幅值电信201');xlabel('以pi为单位的频率');ylabel('绝对值Y');

axis([-1,1,0,60]);

subplot(2,2,4);plot(w/pi,angle(Y)/pi);

title('相位电信201');xlabel('以pi为单位的频率');ylabel('绝对值Y');

axis([-1,1,-1,1]);

b=1;

a=[1,-0.8];

n=0:100;

x=cos(0.05*pi*n);

y=filter(b,a,x);

subplot(2,1,1); stem(n,x)

xlabel('n');ylabel('x(n)');title('输入序列电信1201');

subplot(2,1,2); stem(n,y)

xlabel('n');ylabel('y(n)');title('输出序列电信1201');

实验三

例:已知信号由15Hz幅值0.5的正弦信号和40Hz幅值2的正弦信号组成,数据采样频率fs=100;

N=128;

n=0:N-1;

t=n/fs;

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y=fft(x,N);

f=(0:length(y)-1)'*fs/length(y);

mag=abs(y);

stem(f,mag);

title('N=128点电信1201')

已知带有测量噪声信号f1=50Hz,f2=1为均值为零、方差为1的随机信号,采样频率为1000Hz,t=0:0.001:0.6;

x=sin(2*pi*50*t)+sin(2*pi*120*t);

y=x+2*randn(1,length(t));

Y=fft(y,512);

P=Y.*conj(Y)/512; %求功率

f=1000*(0:255)/512;

subplot(2,1,1);

plot(y);

title('电信1201')

subplot(2,1,2);

plot(f,P(1:256));

title('电信1201')

对信号进行DFT,对其结果进行IDFT,并将IDFT的结果和原信号进行比较。

fs=100; N=128; n=0:N-1; t=n/fs;

x=sin(2*pi*40*t)+sin(2*pi*15*t);

subplot(2,2,1)

plot(t,x)

title('original signal 电信1201')

y=fft(x,N);

mag=abs(y);

f=(0:length(y)-1)'*fs/length(y);

subplot(2,2,2)

plot(f,mag)

title('FFT to original signal 电信1201')

xifft=ifft(y);

magx=real(xifft);

ti=[0:length(xifft)-1]/fs;

subplot(2,2,3)

plot(ti,magx);

title('signal from IFFT 电信1201')

yif=fft(xifft,N);

mag=abs(yif);

subplot(2,2,4)

plot(f,mag)

title('FFT to signal from IFFT 电信201')

用函数conv和FFT计算同一序列的卷积,比较其计算时间。

L=5000; N=L*2-1; n=1:L;

x1=0.5*n; x2=2*n;

t0=clock; yc=conv(x1,x2);

conv_time=etime(clock,t0)

t0=clock;

yf=ifft(fft(x1,N).*fft(x2,N));

fft_time=etime(clock,t0)

实验四

例:用冲激响应不变法设计Butterworth低通数字滤波器,通带波纹小于1dB,阻带在

wp=0.2*pi; ws=0.3*pi;

rp=1;rs=15;ts=0.01;Nn=128;

Wp=wp/ts; Ws=ws/ts;

[N,Wn]=buttord(Wp,Ws,rp,rs,'s');

[z,p,k]=buttap(N);

[Bap,Aap]=zp2tf(z,p,k);

[b,a]=lp2lp(Bap,Aap,Wn);

[bz,az]=impinvar(b,a,1/ts);

freqz(bz,az,Nn,1/ts)

title('电信1201')

设计一个Butterworth高通数字滤波器,通带边界频率为300Hz,阻带边界频率为200Hz,通带波纹小于1dB,阻带衰减大于20dB,采样频率为1000Hz。

fs=1000;

wp=300/(fs/2); ws=200/(fs/2);

rp=1; rs=15; Nn=128;

[N,Wn]=buttord(wp,ws,rp,rs);

[b,a]=butter(N,Wn,'high');

freqz(b,a,Nn,fs)

title('电信1201')

设计一个24阶FIR带通滤波器,通带频率为

wp=[0.35,0.65];

N=24;

b=fir1(2*N,wp);

freqz(b,1,512)

title('电信1201')

f=0:0.002:1;

m(1:201)=1;m(202:301)=0;

m(302:351)=0.5;m(352:401)=0;m(402:501)=1;

hold

plot(f,m,'r:')

b=fir2(64,f,m);

[h,f1]=freqz(b);

f1=f1./pi;

plot(f1,abs(h))

title('N=64 电信1201')

f=0:0.002:1;

m(1:201)=1;m(202:301)=0;

m(302:351)=0.5;m(352:401)=0;m(402:501)=1;

hold

plot(f,m,'r:')

b=fir2(256,f,m);

[h,f1]=freqz(b);

f1=f1./pi;

plot(f1,abs(h))

title('N=256 电信1201')

针对一个含有5Hz、15Hz和30Hz的混和正弦波信号,设计一个FIR带通滤波器。fc1=10; fc2=20; fs=100;

[n,Wn,beta,ftype]=kaiserord([7 13 17 23],[0 1 0], [0.01 0.01

0.01],100);

w1=2*fc1/fs; w2=2*fc2/fs;

window=kaiser(n+1,beta); %使用kaiser窗函数b=fir1(n,[w1 w2],window); %使用标准频率响应的freqz(b,1,512); %数字滤波器频率响应t = (0:100)/fs;

s = sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30);

sf = filter(b,1,s); %对信号s进行滤波figure

subplot(2,1,1); plot(t,s); title(' 电信1201') subplot(2,1,2); plot(t,sf)

title(' 电信1201')

一个简单的Matlab_GUI编程实例

Matlab GUI编程教程(适用于初学者) 1.首先我们新建一个GUI文件:如下图所示; 选择Blank GUI(Default) 2.进入GUI开发环境以后添加两个编辑文本框,6个静态文本框,和一个按钮,布置如下

图所示; 布置好各控件以后,我们就可以来为这些控件编写程序来实现两数相加的功能了。3.我们先为数据1文本框添加代码; 点击上图所示红色方框,选择edit1_Callback,光标便立刻移到下面这段代码的位置。 1. 2. 3.function edit1_Callback(hObject, eventdata, handles) 4.% hObject handle to edit1 (see GCBO) 5.% eventdata reserved - to be defined in a future version of MATLAB

6.% handles structure with handles and user data (see GUIDATA) 7.% Hints: get(hObject,'String') returns contents of edit1 as text 8.% str2double(get(hObject,'String')) returns contents of edit1 as a double 复制代码 然后在上面这段代码的下面插入如下代码: 1. 2.%以字符串的形式来存储数据文本框1的内容. 如果字符串不是数字,则现实空白内容input = str2num(get(hObject,'String')); %检查输入是否为空. 如果为空,则默认显示为0if (isempty(input)) set(hObject,'String','0')endguidata(hObject, handles); 复制代码 这段代码使得输入被严格限制,我们不能试图输入一个非数字。 4.为edit2_Callback添加同样一段代码 5 现在我们为计算按钮添加代码来实现把数据1和数据2相加的目的。 用3中同样的方法在m文件中找到pushbutton1_Callback代码段 如下; 1.function pushbutton1_Callback(hObject, eventdata, handles) 2.% hObject handle to pushbutton1 (see GCBO) 3.% eventdata reserved - to be defined in a future version of MATLAB 4.% handles structure with handles and user data (see GUIDATA) 复制代码

数字信号处理Matlab实现实例(推荐给学生)

数字信号处理Matlab 实现实例 第1章离散时间信号与系统 例1-1 用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB程序如下: a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 图1.1给出了卷积结果的图形,求得的结果存放在数组c中为:{-2 -4 1 3 1 5 1 -3}。 例1-2 用MATLAB计算差分方程 当输入序列为时的输出结果。 解 MATLAB程序如下: N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)];

k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n');ylabel('幅度') 图 1.2 给出了该差分方程的前41个样点的输出,即该系统的单位脉冲响应。 例1-3 用MATLAB 计算例1-2差分方程 所对应的系统函数的DTFT 。 解 例1-2差分方程所对应的系统函数为: 123 123 0.80.440.360.02()10.70.450.6z z z H z z z z -------++= +-- 其DTFT 为 23230.80.440.360.02()10.70.450.6j j j j j j j e e e H e e e e ωωωω ωωω--------++= +-- 用MATLAB 计算的程序如下: k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('实部') xlabel('\omega/\pi');ylabel('幅度')

matlab源代码实例

1.硬币模拟试验 源代码: clear; clc; head_count=0; p1_hist= [0]; p2_hist= [0]; n = 1000; p1 = 0.3; p2=0.03; head = figure(1); rand('seed',sum(100*clock)); fori = 1:n tmp = rand(1); if(tmp<= p1) head_count = head_count + 1; end p1_hist (i) = head_count /i; end figure(head); subplot(2,1,1); plot(p1_hist); grid on; hold on; xlabel('重复试验次数'); ylabel('正面向上的比率'); title('p=0.3试验次数N与正面向上比率的函数图'); head_count=0; fori = 1:n tmp = rand(1); if(tmp<= p2) head_count = head_count + 1; end p2_hist (i) = head_count /i; end figure(head); subplot(2,1,2); plot(p2_hist); grid on; hold on; xlabel('重复试验次数'); ylabel('正面向上的比率'); title('p=0.03试验次数N与正面向上比率的函数图'); 实验结果:

2.不同次数的随机试验均值方差比较 源代码: clear ; clc; close; rand('seed',sum(100*clock)); Titles = ['n=5时' 'n=20时' 'n=25时' 'n=50时' 'n=100时']; Titlestr = cellstr(Titles); X_n_bar=[0]; %the samples of the X_n_bar X_n=[0]; %the samples of X_n N=[5,10,25,50,100]; j=1; num_X_n = 100; num_X_n_bar = 100; h_X_n_bar = figure(1);

MATLAB简单程序大全

MATLAB简单程序大全 求特征值特征向量 A=[2 3 4;1 5 9;8 5 2] det(A) A' rank(A) inv(A) rref(A) eig(A)%求特征值和特征向量 卫星运行问题 h=200,H=51000,R=6378; a=(h+H+2*R)/2; c=(H-h)/2; b=(a^2-c^2)^(1/2); e=c/a; f=sqrt(1-exp(2).*cos(t)^2); l=int(f,t,0,pi/2) L=4*a.*l 动态玫瑰线 n=3;N=10000; theta=2*pi*(0:N)/N; r=cos(n*theta); x=r.*cos(theta); y=r.*sin(theta); comet(x,y) 二重积分 syms x y f=x^2*sin(y); int(int(f,x,0,1),y,0,pi) ezmesh(f,[0,1,0,pi]) 函数画图 syms x;f=exp(-0.2*x)*sin(0.5*x); ezplot(f,[0,8*pi])

玫瑰线 theta=0:0.01:2*pi; r=cos(3*theta); polar(theta,r,'r') 求x^2+y^2=1和x^2+z^2=1所围成的体积 syms x y z R r=1; Z=sqrt(1-x^2); y0=Z; V=8*int(int(Z,y,0,y0),x,0,1) 求导数及图像 f='1/(5+4*cos(x))'; subplot(1,2,1);ezplot(f) f1=diff(f) subplot(1,2,2);ezplot(f1) 绕x轴旋转 t=(0:20)*pi/10; r=exp(-.2*t).*sin(.5*t); theta=t; x=t'*ones(size(t)); y=r'*cos(theta); z=r'*sin(theta); mesh(x,y,z) colormap([0 0 0]) 某年是否闰年 year=input('input year:='); n1=year/4; n2=year/100; n3=year/400; if n1==fix(n1)&n2~=fix(n2) disp('是闰年') elseif n1==fix(n1)&n3==fix(n3) disp('是闰年') else

matlab经典编程例题

以下各题均要求编程实现,并将程序贴在题目下方。 1.从键盘输入任意个正整数,以0结束,输出那些正整数中的素数。 clc;clear; zzs(1)=input('请输入正整数:');k=1; n=0;%素数个数 while zzs(k)~=0 flag=0;%是否是素数,是则为1 for yz=2:sqrt(zzs(k))%因子从2至此数平方根 if mod(zzs(k),yz)==0 flag=1;break;%非素数跳出循环 end end if flag==0&zzs(k)>1%忽略0和1的素数 n=n+1;sus(n)=zzs(k); end k=k+1; zzs(k)=input('请输入正整数:'); end disp(['你共输入了' num2str(k-1) '个正整数。它们是:']) disp(zzs(1:k-1))%不显示最后一个数0 if n==0 disp('这些数中没有素数!')%无素数时显示 else disp('其中的素数是:') disp(sus) end 2.若某数等于其所有因子(不含这个数本身)的和,则称其为完全数。编程求10000以内所有的完全数。 clc;clear;

wq=[];%完全数赋空数组 for ii=2:10000 yz=[];%ii的因子赋空数组 for jj=2:ii/2 %从2到ii/2考察是否为ii的因子 if mod(ii,jj)==0 yz=[yz jj];%因子数组扩展,加上jj end end if ii==sum(yz)+1 wq=[wq ii];%完全数数组扩展,加上ii end end disp(['10000以内的完全数为:' num2str(wq)])%输出 3.下列这组数据是美国1900—2000年人口的近似值(单位:百万)。 (1)若. 2c + = y+ 与试编写程序计算出上式中的a、b、c; 的经验公式为 t at bt y (2)若.bt 的经验公式为 y= 与试编写程序计算出上式中的a、b; y ae t (3)在一个坐标系下,画出数表中的散点图(红色五角星),c + =2中 ax bx y+拟合曲线图(蓝色实心线),以及.bt y=(黑色点划线)。 ae (4)图形标注要求:无网格线,横标注“时间t”,纵标注“人口数(百万)”,图形标题“美国1900—2000年的人口数据”。 (5)程序中要有注释,将你的程序和作好的图粘贴到这里。 clf;clc;clear %清除图形窗、屏幕、工作空间 t=1900:10:2000; y=[76 92 106 123 132 151 179 203 227 250 281]; p1=polyfit(t,y,2);%二次多项式拟合

三个遗传算法matlab程序实例

遗传算法程序(一): 说明: fga.m 为遗传算法的主程序; 采用二进制Gray编码,采用基于轮盘赌法的非线性排名选择, 均匀交叉,变异操作,而且还引入了倒位操作! function [BestPop,Trace]=fga(FUN,LB,UB,eranum,popsize,pCross,pMutation,pInversion,options) % [BestPop,Trace]=fmaxga(FUN,LB,UB,eranum,popsize,pcross,pmutation) % Finds a maximum of a function of several variables. % fmaxga solves problems of the form: % max F(X) subject to: LB <= X <= UB % BestPop - 最优的群体即为最优的染色体群 % Trace - 最佳染色体所对应的目标函数值 % FUN - 目标函数 % LB - 自变量下限 % UB - 自变量上限 % eranum - 种群的代数,取100--1000(默认200) % popsize - 每一代种群的规模;此可取50--200(默认100) % pcross - 交叉概率,一般取0.5--0.85之间较好(默认0.8) % pmutation - 初始变异概率,一般取0.05-0.2之间较好(默认0.1) % pInversion - 倒位概率,一般取0.05-0.3之间较好(默认0.2) % options - 1*2矩阵,options(1)=0二进制编码(默认0),option(1)~=0十进制编 %码,option(2)设定求解精度(默认1e-4) % % ------------------------------------------------------------------------ T1=clock; if nargin<3, error('FMAXGA requires at least three input arguments'); end if nargin==3, eranum=200;popsize=100;pCross=0.8;pMutation=0.1;pInversion=0.15;options=[0 1e-4];end if nargin==4, popsize=100;pCross=0.8;pMutation=0.1;pInversion=0.15;options=[0 1e-4];end if nargin==5, pCross=0.8;pMutation=0.1;pInversion=0.15;options=[0 1e-4];end if nargin==6, pMutation=0.1;pInversion=0.15;options=[0 1e-4];end if nargin==7, pInversion=0.15;options=[0 1e-4];end if find((LB-UB)>0) error('数据输入错误,请重新输入(LB

MATLAB程序设计及经典例题解析3

MATLAB程序设计 用MATLAB语言编写的程序,称为M文件。M文件可以根据调用方式的不同分为两类:命令文件(Script File)和函数文件(Function File)。 例3-1 分别建立命令文件和函数文件,将华氏温度f转换为摄氏温度c。 程序1:首先建立命令文件并以文件名f2c.m存盘。 clear; %清除工作空间中的变量 f=input('Input Fahrenheit temperature:'); c=5*(f-32)/9 然后在MATLAB的命令窗口中输入f2c,将会执行该命令文件,执行情况为: Input Fahrenheit temperature:73 c =22.7778 程序2:首先建立函数文件f2c.m。 function c=f2c(f) c=5*(f-32)/9 然后在MATLAB的命令窗口调用该函数文件。 clear; y=input('Input Fahrenheit temperature:'); x=f2c(y) 输出情况为: Input Fahrenheit temperature:70 c =21.1111 x =21.1111 3.1.2 M文件的建立与打开 M文件是一个文本文件,它可以用任何编辑程序来建立和编辑,而一般常用且最为方便的是使用MATLAB提供的文本编辑器。

1.建立新的M文件 为建立新的M文件,启动MATLAB文本编辑器有3种方法: (1) 菜单操作。从MATLAB主窗口的File菜单中选择New菜单项,再选择M-file命令,屏幕上将出现MATLAB 文本编辑器窗口。 (2) 命令操作。在MATLAB命令窗口输入命令edit,启动MATLAB文本编辑器后,输入M文件的内容并存盘。 (3) 命令按钮操作。单击MATLAB主窗口工具栏上的New M-File命令按钮,启动MATLAB文本编辑器后,输入M文件的内容并存盘。 2.打开已有的M文件 打开已有的M文件,也有3种方法: (1) 菜单操作。从MATLAB主窗口的File菜单中选择Open命令,则屏幕出现Open对话框,在Open对话框中选中所需打开的M文件。在文档窗口可以对打开的M文件进行编辑修改,编辑完成后,将M文件存盘。 (2) 命令操作。在MATLAB命令窗口输入命令:edit 文件名,则打开指定的M文件。 (3) 命令按钮操作。单击MATLAB主窗口工具栏上的Open File命令按钮,再从弹出的对话框中选择所需打开的M文件。 3.2 程序控制结构 3.2.1 顺序结构 1.数据的输入 从键盘输入数据,则可以使用input函数来进行,该函数的调用格式为: A=input(提示信息,选项); 其中提示信息为一个字符串,用于提示用户输入什么样的数据。 如果在input函数调用时采用's'选项,则允许用户输入一个字符串。例如,想输入一个人的姓名,可采用命令: xm=input('What''s your name?','s'); 2.数据的输出 MATLAB提供的命令窗口输出函数主要有disp函数,其调用格式为

图论算法及matlab程序的三个案例

图论实验三个案例 单源最短路径问题 Dijkstra 算法 Dijkstra 算法是解单源最短路径问题的一个贪心算法。其基本思想是,设置一个顶点集合S 并不断地作贪心选择来扩充这个集合。一个顶点属于集合S 当且仅当从源到该顶点的最短路径长度已知。设v 是图中的一个顶点,记()l v 为顶点v 到源点v 1的最短距离,,i j v v V ?∈,若 (,)i j v v E ?,记i v 到j v 的权ij w =∞。 Dijkstra 算法: ① 1{}S v =,1()0l v =;1{}v V v ??-,()l v =∞,1i =,1{}S V v =-; ② S φ=,停止,否则转③; ③ ()min{(),(,)} j l v l v d v v =, j v S ∈,v S ?∈; ④ 存在1 i v +,使 1()min{()} i l v l v +=,v S ∈; ⑤ 1{} i S S v +=U , 1{} i S S v +=-,1i i =+,转②; 实际上,Dijkstra 算法也是最优化原理的应用:如果121n n v v v v -L 是从1v 到 n v 的最 短路径,则 121 n v v v -L 也必然是从1v 到 1 n v -的最优路径。 在下面的MATLAB 实现代码中,我们用到了距离矩阵,矩阵第i 行第j 行元素表 示顶点i v 到j v 的权ij w ,若i v 到j v 无边,则realmax ij w =,其中realmax 是MATLAB 常量,表示最大的实数+308)。 function re=Dijkstra(ma) %用Dijkstra 算法求单源最短路径 %输入参量ma 是距离矩阵 %输出参量是一个三行n 列矩阵,每列表示顶点号及顶点到源的最短距离和前顶点 n=size(ma,1);%得到距离矩阵的维数 s=ones(1,n);s(1)=0;%标记集合S 和S 的补 r=zeros(3,n);r(1,:)=1:n;r(2,2:end)=realmax;%初始化 for i=2:n;%控制循环次数 mm=realmax; for j=find(s==0);%集合S 中的顶点 for k=find(s==1);%集合S 补中的顶点

Matlab简单实例学习

Matlab 程序代码 绘制y = 10e-1.5t sin( 7.75t ) 的函数图象 7.75 fv clear; t=0:0.02:10; f1=10/sqrt(7.75).*exp(-1.5*t); f2=sin(sqrt(7.75).*t); y=f1.*f2; plot(t,y,'-k',t,y,'ok'); xlabel('t');ylabel('y(t) ');title('函数图像') axis([-2 10 -0.5 2]) 拉氏变换 clear; clc; syms s t fs1 fs2 fs3 ft1 ft2 ft3; L=1,C=0.1,R=[1.5 3 5]; h1=1/(L*C*s^2+R(1)*C*s+1); h2=1/(L*C*s^2+R(2)*C*s+1); h3=1/(L*C*s^2+R(3)*C*s+1); fs1=h1*(1/s); fs2=h2*(1/s); fs3=h3*(1/s);

ft1=ilaplace(fs1,s,t); ft2=ilaplace(fs2,s,t); ft3=ilaplace(fs3,s,t); ezplot(t,ft1); hold on; ezplot(t,ft2); hold on; ezplot(t,ft3); 信号编码 对[1 1 0 1 1 1 0 1 0 0 1]进行编码。 clear; clc; c=[1 1 0 1 1 1 0 1 0 0 1] for i=1:length(c) if i==1 d1(i)=0;d2(i)=0; elseif i==2 d1(i)=c(i-1);d2(i)=c(i-1); elseif i==3 d1(i)=mod(c(i-1)+c(i-2),2); d2(i)=c(i-1); else d1(i)=mod(c(i-1)+c(i-2),2); d2(i)=mod(c(i-1)+c(i-3),2); end

矩量法matlab程序设计实例

矩量法m atla b程序设计实例: Ha llen 方程求对称振子天线 一、条件与计算目标 已知: 对称振子天线长为L,半径为a ,且天线长度与波长得关系为,,设,半径a=0、0000001,因此波数为。 目标: 用H all en 方程算出半波振子、全波振子以及不同值得对应参数值。 求:(1)电流分布 (2)E 面方向图 (二维),H 面方向图(二维),半波振子空间方向性图(三维) 二、对称振子放置图 图1 半波振子得电流 分布 半波振子天线平行于z 轴放置,在x轴与y轴上得分量都为零,坐标选取方式有两种形式,一般选取图1得空间放置方 式。图1给出了天线得电流分布情况,由图可知,当天线很细时,电流分布近似正弦分布。 三、Ha llen 方程 得解题思路 ()()()()2 1 ' ' ' ' 12,cos sin sin 'z z i z z z z i z k z G z z dz c kz c kz E k z z dz j ωμ'++=-?? 对于中心馈电得偶极子,Hallen 方程为 ()22'1222 ('),'cos sin sin ,2L L i L L V i z G z z dz c kz c kz k z z j η + -- ++= <<+? 脉冲函数展开与点选配,得到 ()1121 ,''cos sin sin ,1,2,,2n n N z i n m m m m z n V I G z z dz c kz c kz k z m N j η +''=++= =???∑? 上式可以写成 矩阵形式为 四、结果与分析 (1)电流分布

matlab程序设计实例

MATLAB 程序设计方法及若干程序实例 樊双喜 (河南大学数学与 信息科学学院开封475004) 摘要本文通过对 MATLAB 程序设计中的若干典型问题做简要的分析和总结,并在此基础上着重讨论了有关算法设计、程序的调试与测试、算法与程序的优化以及循环控制等方面的问题.还通过对一些程序实例做具体解析,来方便读者进行编程训练并掌握一些有关MATLAB 程序设计方面的基本概念、基本方法以及某些问题的处理技巧等.此外,在文章的最后还给出了几个常用数学方法的算法程序, 供读者参考使用.希望能对初学者进行 MATLAB 编程训练提供一些可供参考的材料,并起到一定的指导和激励作用,进而为MATLAB 编程入门打下好的基础. 关键字算法设计;程序调试与测试;程序优化;循环控制 1 算法与程序 1.1 算法与程序的关系算法被称为程序的灵魂,因此在介绍程序之前应先了 解什么是算法.所谓算 法就是对特定问题求解步骤的一种描述.对于一个较复杂的计算或是数据处理的问题,通常是先设计出在理论上可行的算法,即程序的操作步骤,然后再按照算法逐步翻译成相应的程序语言,即计算机可识别的语言. 所谓程序设计,就是使用在计算机上可执行的程序代码来有效的描述用于解决特定问题算法的过程.简单来说,程序就是指令的集合.结构化程序设计由于采用了模块分化与功能分解,自顶向下,即分而治之的方法,因而可将一个较复杂的问题分解为若干子问题,逐步求精.算法是操作的过程,而程序结构和程序流程则是算法的具体体现. 1.2MATLAB 语言的特点 MATLAB 语言简洁紧凑,使用方便灵活,库函数极其丰富,其语法规则与科技人员的思维和书写习惯相近,便于操作.MATLAB 程序书写形式自由,利用其丰富

matlab程序设计例题及答案

1.编写程序:计算1/3+2/5+3/7+……+10/21 法一: s=0; for i=1:10 s=s+i/(2*i+1); end s s = 4.4096 法二: sum((1:10)./(3:2:21)) ans = 4.4096 2.编写程序:计算1~100中即能被3整除,又能被7整除的所有数之和。 s=0; for i=1:100 if mod(i,3)==0&&mod(i,7)==0 s=s+i; end,end s s = 210 3.画出y=n!的图(1<=n<=10),阶乘的函数自己编写,禁用MATLAB自带的阶乘函数。 x=1:10; for i=1:10 try y(i)=y(i-1)*i; catch y(i)=1; end,end plot(x,y)

12345678910 0.511.522.533.54x 10 6 4.一个数恰好等于它的因子之和,这个数就称为完数。例如,6的因子为1,2,3,而6=1+2+3,因此6就是一个完数。编程找出2000以内的所有完数。 g=[]; for n=2:2000 s=0; for r=1:n-1 if mod(n,r)==0 s=s+r; end end if s==n g=[g n]; end end g g =6 28 496

5.编写一个函数,模拟numel函数的功能,函数中调用size函数。 function y=numelnumel(x) m=size(x); y=m(1)*m(2); numelnumel([1 2 3;4 5 6]) ans = 6 6. 编写一个函数,模拟length函数的功能,函数中调用size函数。 function y=lengthlength(x) m=size(x); y=max(m(1),m(2)); lengthlength([1 2 3;4 5 6]) ans = 3 7.求矩阵rand(5)的所有元素和及各行平均值,各列平均值。 s=rand(5); sum=sum(sum(s)) mean2=mean(s,2) mean1=mean(s) sum = 13.8469

多目标优化实例和matlab程序

NSGA-II 算法实例 目前的多目标优化算法有很多,Kalyanmoy Deb 的带精英策略的快速非支配排序遗传算法(NSGA-II)无疑是其中应用最为广泛也是最为成功的一种。本文用的算法是MATLAB 自带的函数gamultiobj ,该函数是基于NSGA-II 改进的一种多目标优化算法。 一、数值例子 多目标优化问题 42422 11211122124224212212112 12min (,)10min (,)55..55 f x x x x x x x x x f x x x x x x x x x s t x =-++-=-++-≤≤??-≤≤?二、Matlab 文件 1.适应值函数m 文件: function y=f(x) y(1)=x(1)^4-10*x(1)^2+x(1)*x(2)+x(2)^4-x(1)^2*x(2)^2; y(2)=x(2)^4-x(1)^2*x(2)^2+x(1)^4+x(1)*x(2);2.调用gamultiobj 函数,及参数设置: clear clc fitnessfcn=@f; %适应度函数句柄nvars=2; %变量个数lb=[-5,-5]; %下限ub=[5,5]; %上限A=[];b=[];%线性不等式约束 Aeq=[];beq=[];%线性等式约束 options=gaoptimset('paretoFraction',0.3,'populationsize',100,'generations',200,'stallGenLimit',200,'TolFun',1e-100,'PlotFcns',@gaplotpareto); %最优个体系数paretoFraction 为0.3;种群大小populationsize 为100,最大进化代数generations 为200, %停止代数stallGenLimit 为200,适应度函数偏差TolFun 设为1e-100,函数gaplotpareto :绘制Pareto 前端 [x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)

图论算法及matlab程序的三个案例

图论实验三个案例 单源最短路径问题 1.1 Dijkstra 算法 Dijkstra 算法是解单源最短路径问题的一个贪心算法。其基本思想是,设置一个顶点集合S 并不断地作贪心选择来扩充这个集合。一个顶点属于集合S 当且仅当从源到该顶点的最短路径长度已知。设v 是图中的一个顶点,记()l v 为顶点 v 到源点v 1的最短距离, ,i j v v V ?∈,若 (,)i j v v E ?,记i v 到 j v 的权 ij w =∞ 。 Dijkstra 算法: ① 1{}S v =,1()0l v =;1{}v V v ??-,()l v =∞,1i =,1{}S V v =-; ② S φ=,停止,否则转③; ③ ()min{(),(,)} j l v l v d v v =, j v S ∈,v S ?∈; ④ 存在1i v +,使1()min{()}i l v l v +=,v S ∈; ⑤ 1{}i S S v += ,1{}i S S v +=-,1i i =+,转②; 实际上,Dijkstra 算法也是最优化原理的应用:如果121n n v v v v - 是从1v 到n v 的最短路径,则121n v v v - 也必然是从1v 到1n v -的最优路径。 在下面的MATLAB 实现代码中,我们用到了距离矩阵,矩阵第i 行第j 行元素表示顶点i v 到 j v 的权 ij w ,若i v 到 j v 无边,则 realmax ij w =,其中realmax 是 MATLAB 常量,表示最大的实数(1.7977e+308)。 function re=Dijkstra(ma)

Matlab简单实例学习

Matlab 程序代码 绘制 1.5 10sin(7.75)7.75 t y e t -= 的函数图象。 fv clear; t=0:0.02:10; f1=10/sqrt(7.75).*exp(-1.5*t); f2=sin(sqrt(7.75).*t); y=f1.*f2; plot(t,y,'-k',t,y,'ok'); xlabel('t');ylabel('y(t) ');title('函数图像') axis([-2 10 -0.5 2]) 拉氏变换 clear; clc; syms s t fs1 fs2 fs3 ft1 ft2 ft3; L=1,C=0.1,R=[1.5 3 5]; h1=1/(L*C*s^2+R(1)*C*s+1); h2=1/(L*C*s^2+R(2)*C*s+1);

h3=1/(L*C*s^2+R(3)*C*s+1); fs1=h1*(1/s); fs2=h2*(1/s); fs3=h3*(1/s); ft1=ilaplace(fs1,s,t); ft2=ilaplace(fs2,s,t); ft3=ilaplace(fs3,s,t); ezplot(t,ft1); hold on; ezplot(t,ft2); hold on; ezplot(t,ft3); 信号编码 对[1 1 0 1 1 1 0 1 0 0 1]进行编码。 clear; clc; c=[1 1 0 1 1 1 0 1 0 0 1] for i=1:length(c) if i==1 d1(i)=0;d2(i)=0; elseif i==2 d1(i)=c(i-1);d2(i)=c(i-1); elseif i==3

遗传算法的MATLAB程序实例讲解学习

遗传算法的M A T L A B 程序实例

遗传算法的程序实例 如求下列函数的最大值 f(x)=10*sin(5x)+7*cos(4x) x∈[0,10] 一、初始化(编码) initpop.m函数的功能是实现群体的初始化,popsize表示群体的大小,chromlength表示染色体的长度(二值数的长度), 长度大小取决于变量的二进制编码的长度(在本例中取10位)。 代码: %Name: initpop.m %初始化 function pop=initpop(popsize,chromlength) pop=round(rand(popsize,chromlength)); % rand随机产生每个单元为 {0,1} 行数为popsize,列数为chromlength的矩阵, % roud对矩阵的每个单元进行圆整。这样产生的初始种群。 二、计算目标函数值 1、将二进制数转化为十进制数(1) 代码: %Name: decodebinary.m %产生 [2^n 2^(n-1) ... 1] 的行向量,然后求和,将二进制转化为十进制 function pop2=decodebinary(pop) [px,py]=size(pop); %求pop行和例数 for i=1:py pop1(:,i)=2.^(py-1).*pop(:,i); py=py-1; end pop2=sum(pop1,2); %求pop1的每行之和 2、将二进制编码转化为十进制数(2) decodechrom.m函数的功能是将染色体(或二进制编码)转换为十进制,参数spoint表示待解码的二进制串的起始位置。(对于多个变量而言,如有两个变量,采用20为表示,每个变量10为,则第一个变量从1开始,另一个变量从11开始。本例为1),参数1ength表示所截取的长度(本例为10)。 代码: %Name: decodechrom.m %将二进制编码转换成十进制 function pop2=decodechrom(pop,spoint,length) pop1=pop(:,spoint:spoint+length-1); pop2=decodebinary(pop1); 3、计算目标函数值 calobjvalue.m函数的功能是实现目标函数的计算,其公式采用本文示例仿真,可根据不同优化问题予以修改。

matlab30个案例分析案例6代码

Draw %function J=draw(individual) load best zbest individual=zbest; %函数功能:画出最优粒子对应的各种图形 %individual输入粒子 %fitness输出适应度值 w11=reshape(individual(1:6),3,2); w12=reshape(individual(7:12),3,2); w13=reshape(individual(13:18),3,2); w21=individual(19:27); w22=individual(28:36); w23=individual(37:45); rate1=0.006;rate2=0.001;%学习率 k=0.3;K=3; y_1=zeros(3,1);y_2=y_1;y_3=y_2;%输出值 u_1=zeros(3,1);u_2=u_1;u_3=u_2;%控制率 h1i=zeros(3,1);h1i_1=h1i;%第一个控制量 h2i=zeros(3,1);h2i_1=h2i;%第二个控制量 h3i=zeros(3,1);h3i_1=h3i;%第三个空置量 x1i=zeros(3,1);x2i=x1i;x3i=x2i;x1i_1=x1i;x2i_1=x2i;x3i_1=x3i;%隐含层输出 %权值初始化 k0=0.03; %值限定 ynmax=1;ynmin=-1;%系统输出值限定 xpmax=1;xpmin=-1;%P节点输出限定 qimax=1;qimin=-1;%I节点输出限定 qdmax=1;qdmin=-1;%D节点输出限定 uhmax=1;uhmin=-1;%输出结果限定 for k=1:1:200 %--------------------------------网络前向计算-------------------------- %系统输出 y1(k)=(0.4*y_1(1)+u_1(1)/(1+u_1(1)^2)+0.2*u_1(1)^3+0.5*u_1(2))+0.3*y_1(2); y2(k)=(0.2*y_1(2)+u_1(2)/(1+u_1(2)^2)+0.4*u_1(2)^3+0.2*u_1(1))+0.3*y_1(3); y3(k)=(0.3*y_1(3)+u_1(3)/(1+u_1(3)^2)+0.4*u_1(3)^3+0.4*u_1(2))+0.3*y_1(1);

Matlab实现Zoutendijk编程例子

22 212312132312min f (x)x 2x 3x x x 2x x x x 4x 6x =+++-+-- . 123123x 2x x 4 x ,x ,x 0++≤≥ 取初始点(1)T x (0,0,0)=,通过Matlab 编程实现求解过程。 公用函数如下: 1、function [f,x]=func %设置目标函数 syms x1 x2 x3; f=x1^2+2*x2^2+3*x3^2+x1*x2-2*x1*x3+x2*x3-4*x1-6*x2; x=[x1,x2,x3]; end 2、function f_val=fval(x0) %求目标函数值 x0=transpose(x0); [f,x]=func; f_val=subs(f,x,x0); end 3、function s=diff_val(x0) %求目标函数梯度 [f,x]=func; grad=jacobian(f,x); s=subs(grad,x,x0); end 4、function h=fmin(x0,d0,vmax) %求函数最小值 [f,x]=func; syms h ; a=x0+h*d0; f_val=inline(subs(f,x,a)); if vmax==inf min_h=fminbnd(f_val,0,10000); else min_h=fminbnd(f_val,0,vmax); end h=min_h; end Zoutendijk 方法主函数 function [X0,f_val]=zoutendijk(A,b,x0,Aeq,beq) %自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数

matlab程序设计例题及答案

matlab程序设计例题及答案 1.编写程序:计算1/3+2/5+3/7+……+10/21 法一: s=0; for i=1:10 s=s+i/(2*i+1); end ss = 法二: sum((1:10)./(3:2:21)) ans = 2.编写程序:计算1~100中即能被3整除,又能被7整除的所有数之和。 s=0; for i=1:100 if mod(i,3)==0&&mod(i,7)==0 s=s+i; end,end ss = 210 3.画出y=n!的图,阶乘的函数自己编写,禁用MATLAB 自带的阶乘函数。 x=1:10; for i=1:10 try y(i)=y(i-1)*i; catch y(i)=1; end,end plot(x,y) 10612345678910 4.一个数恰好等于它的因子之和,这个数就称为完数。

例如,6的因子为1,2,3,而6=1+2+3,因此6就是一个完数。编程找出20XX以内的所有完数。 g=; for n=2:20XX s=0; for r=1:n-1 if mod(n,r)==0 s=s+r; end end if s==n g=[g n]; end end g g =6 28 496 5.编写一个函数,模拟numel函数的功能,函数中调用size函数。 function y=numelnumel(x) m=size(x); y=m(1)*m(2); numelnumel([1 2 3;4 5 6]) ans = 6 6. 编写一个函数,模拟length函数的功能,函数中调用size函数。 function y=lengthlength(x) m=size(x); y=max(m(1),m(2)); lengthlength([1 2 3;4 5 6]) ans = 3

相关主题