matlab四合一(dream)

第一题

目的：

1） 学会使用矩阵和向量的思维方式来考虑问题。这是因为Matlab以该种方式组织数据，通

过这种方式，一个复杂运算在选择了合适的函数和数据组织方式基础上，往往可以通过1～2行代码解决。

2） 学会和习惯使用help来学习新的函数用法。需要用到的函数本文在题目后面标出。

题目说明：问题1-7请新建一个脚本，命名为shortProblems.m，7道小题目之间以注释隔开。 B303

1.标量变量生成，请生成如下变量：

1)a?102)b?2.5?10233)c?2?3（ii表示复数虚部）

2.向量生成，请生成如下向量变量：

/34)d?ej2?（j表示复数虚部，用到exp和pi）1)aVec??3.1415926??2.71??8??2)bVec???28???182??3)cVec??54.8...?4.8?5?(5~?5.步长为-0.2)4)dVec?[100100.01...100.99101](1~10确保length正确，用到logsapce)5)eVec?Hello（一个字符串）

3.矩阵变量，请生成如下矩阵变量：

?2...2??(9?9,值为2,用到ones或者zeros)1)aMat?????2??2?0??10...?0?0???(9?9，对角值为[1 2 3 4 5 4 3 2 1],其余值为0，用到diag)2)bMat=?050??00???01??0??111...91??212?92?(10?10,1:100，用到reshape)3)cMat??????1020100???NaNNaNNaNNaN??(3?4,用到nan)4)dMat??NaNNaNNaNNaN????NaNNaNNaNNaN???13?15?5)eMat?????2210?87?6)fMat,一个5?3矩阵，元素值为随机正数，值的范围为-3～3，使用rand，floor或者ceil

4.标量方程，使用题目1中生成的变量计算x，y和z值：

1）x?11?e(?(a?15)/6)g2)y?(a?21b)?(记得吗？h=h1/g，使用sqrt)log(?[(c?d)(c?d)]sin(a?/3))3)z=(?表示复数的实部，c表示c的共轭，使用real，conj，log)cc

5.向量方程计算，使用题目2中生成的变量求解以下方程的值。使用元素对应法则，.*，./，.^。

1)xVec=12?2.52e?cVec2/(2?2.52)2)yVec?(aVecT)2?bVec2(aVecT表转置)3)zVec?log10(1/dVec)(log10表示以10为底的对数，使用log10)

6.矩阵运算，使用题目2和3中生成的变量求解以下方程：

1)xMat?(aVec?bVec)?aMat22)yMat?bVec?aVec(注意结果不同于aVec?bVec)3)zMat?cMat(aMat?bMat)T(cMat表示求cMat的行列式，用到det)

7.一般运算及元素索引

1)按列规则计算cMat的cSum，结果应为一个行向量，使用sum2)按行计算eMat的平均eMean，结果为一个列向量，使用mean3)将eMat的第一行值替换为?111?4)将cMat的第2行～9行和第2列～9列提出，赋给cSub5)已知一个向量为lin=?12...20?,将其中的偶数项变为负值得到lin=?1-23-4将这些值置为0-20?6)使用rand生成一个1?5的向量，记为r，使用find找到其中值小于0.5的，并

shortProblem

clear all; a=10;

b=2.5*10^23; c=2+3*i;

d=exp((j*2*pi)/3);

%%%%%%%%%%%%%%%%%%%%%% aVec=[3.14 15 9 26]; bVec=ceil[2.71 8 28 182]; cVec=length[5 -0.2 5]; dVec=logspace[0 1 101]; eVec='Hello';

%%%%%%%%%%%%%%%%%%%%%%%% n=ones(9) aMat=n*2 bMat=zeros(9)

diag([1 2 3 4 5 4 3 2 1]) cMat=reshape(p,10,10) dMat=nan(3 4)

eMat=[13 -1 5;-22 10 -87] fMat=round(rand(5,3)*6)-3

%%%%%%%%%%%%%%%%%%%%%%%%%%% x=1/(exp(-(a-15)/6)) y=(a^(1/2)+b^(1/21))^pi

z=(log(real[(c+d)*(c-d)]*sin(a*pi/3))/( conj*con(-j))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xVec=(1/((2*pi*2.5^2)^0.5))*exp((-cVec^2)/(2*2.5^2)) yVec=((aVec')^2+b*Vec^2)^0.5 zVec=log10(1/dVec)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xMat=(aVec*bVec)*aMat^2 yMat=bMat*aMat

zMat=det(cMat)*(aMat*bMat)'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sum=cSum(cMat) mean=eMean(eMat)' eMat=[1 1 1;-22 10 -87] cSub=cMat([2 9;9 2]) lin=[1 1 20]

lin=-1*lin[2 2 20] find=(r(1 5)<0.5) find=0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

8.绘制曲线练习

1)新建一个脚本，命名为twoLinePlot.m2)使用figure打开一个绘图窗口3)绘制正弦和余弦一个周期的曲线a.时间向量t为0～2?，自行定义点数以确保采样点足够，曲线才能够平滑b.绘制sin(t)c.使用hold on命令，这会告诉figure不要丢弃之前绘制的曲线，同样可以使用hold off来关闭hold属性d.绘制cos(t)，使用红色，dashed线样式4)添加labelsa.xlabel，ylabelb.使用title给figure一个标题c.使用legend，注明两条曲线各自的名称5)使用xlim和ylim使得曲线看起来更漂亮，xlim设置x轴标度为0～2?，ylim设置y轴标度为-1.4～1.46)运行该脚本结果参考下图：

Sin and Cos SinCos10.5Function value0-0.5-1 0123Time(s)456

twoLinePlot

figure

t1=0:0.01:2*pi; y1=sin(t1); plot(t1,y1); axis([0 6 -1 1])

line([0 6],[0 0],'color','k') line([0 0],[-1 1],'color','k')

xlabel('Time(s)'),ylabel('Function value'); title(['Sin and Cos']) legend('sin','cos',0) hold on;

%%%%%%%%%%%%%%%%%%%%%%%%% t1=0:0.01:2*pi; y2=cos(t1); plot(t1,y2,'r--') axis([0 6 -1 1])

line([0 6],[0 0],'color','k') line([0 0],[-1 1],'color','k')

xlabel('Time(s)'),ylabel('Function value'); title(['Sin and Cos']) legend('sin','cos',0) hold off;

9.变量操作

1)新建一个脚本，命名为calculateGrades.m2)load数据classGrades.mat3)查看载入的数据，输入前5行到命令行。第一列表示学号1～15，其余7列表示7次题目的分数（0～5分），有一些为NaN，表示该同学缺席或者没交作业，我希望你们当中没有4)我们只关心成绩，提取载入数据中的分数（2～8列），赋给变量grades，需要注意的是，请使用end或者size(namesAndGrades,2),这样我们的代码适用于任意维度的矩阵5)计算所有人每次成绩的平均成绩，结果应为1?7的向量，使用mean函数，结果显示在命令行，我们发现含有NaN值的那一列得到结果也为NaN，为了避免这个问题，使用nanmean函数重新计算。将这个结果赋给变量meanGrades，并显示在命令行查看以确保当中没有NaN值。6）设定百分制70分可以评为等级B，那么5分制的话，就是3.5分可以评定为B。a.建立meanMatrix变量，维度跟grades相同，其中每行的值分别是meanGrades中不同列的值，使用repmat来建立。输出到命令行来检查是否正确b.按照公式curvedGrades=3.5（grades/meanMatrix）计算标准化后的成绩，记在变量curvedGrades中，记得要使用./c.使用nanmean计算curvedGrades的平均值查看是否都是3.5d.因为我们的计算方法，有可能一些分数接近5分，那么计算的结果将大于5，采用find函数，将大于5的结果置为57)按列计算curvedGrades的平均结果，使用nanmean，注意维度，得到每个同学7次作业的平均成绩，记为totalGrade。此外，我们希望结果为1～5的正数，因此记得使用ceil。接下来生命一个变量letters，其值为FDCBA。按照totalGrade中的值分别对应等级评定，给出每位同学的最终等级评定结果。使用switch函数，例如成绩为1的，则评为F，成绩为5则为A。最后使用disp函数，输出Grades： letterGrades运行结果应是这样的：

calculateGrades

clear all;

x=xlsread('testdata.xls') ans=x(1:5,:)

grades=x(:,2:end); ans=mean(grades)

meanGrades=nanmean(grades)

meanMatrix=repmat(meanGrades,size(grades,1),1) curvedGrades=3.5*(grades./meanMatrix); ans=nanmean(curvedGrades)

curvedGrades(find(curvedGrades>5))=5; totalGrades=ceil(nanmean(curvedGrades'))' letters='FDCBA';

for k=1:length(totalGrades) switch totalGrades(k) case 1

Grades(k)=letters(1); case 2

Grades(k)=letters(2); case 3

Grades(k)=letters(3); case 4

Grades(k)=letters(4); case 5

Grades(k)=letters(5); end end

disp(['Grades:',Grades]);

10.抛球模型

1）建立一个新脚本，命名问throwBall.m。

2）定义常量：球的初始高度h=1.5m；重力加速度g=9.8m/s2；初始球速为v=4m/s；初始抛球的角度为θ=45度。

3）生成时间向量，0～1之间，1000个点线性分布，使用linspace。

4）我们已知抛球的物理模型：假定x为距离，y为高度，则根据高中物理所学水平位移为

x(t)?vcos(??180)t，垂直位移为y(t)?h?vsin(??180)t?12gt，利用该公式计算x和y。 25）计算何时小球撞倒地面，使用find找到高度变负的地方，那么水平距离就是那时间点的x值，使用disp输出结果：小球在x米处撞倒地面。

6）使用figure打开绘图窗口，绘制y和x的对应曲线，xlabel，ylabel和title，hold on figure，绘制地面位置，使用黑色，虚线。这应该是一个水平曲线， 7）运行代码查看结果。

throwBall

clear all; figure(2) h=1.5; g=9.8; v=4; q=45; q=45

t=linspace(0,1,1000); x=v*cos(q*pi/180)*t; plot(t,x,'r'); hold on;

t=linspace(0,1,1000);

y=h+v*sin(q*pi/180)*t-0.5*g*t.^2; plot(t,y);

t0=min(find(y<=0)); x0=x(t0);

xlabel('时间轴'); ylabel('函数'); title('抛球模型'); legend('x(t)','y(t)');

disp(['小球在',num2str(x0),'处撞到地面']); hold on;

plot([0 max(x)],[0 0],'k--'); hold off;

11.加密算法

1）新建一个脚本命名为encrypt.m。

2）定义原始密码为“This is my top secret message!“。 3）让我们随机对原始密码中的字符进行随机排序，以生成加密文本，使用randperm和length。 4）对原始密码利用加密文本进行加密，加密后的文本命名为encoded。

5）生成一个两列的矩阵，第一列数据为加密文本的序号，第二列为正常的序列从1开始到文本长度，使用length。

6）使用sortrows函数按照第一列对第二列进行排序，以得到解密的序号

7）将第二列序号取出，对于加密后的文本encoded按照解密后的序号进行值的索引则得到解密文本

8）输入原始，加密和解密文本到命令行

9）使用strcmp函数对原始密码和解密文本进行比较

10）解密正确值为1，错误为0，使用disp函数输出结果到命令行模式。 结果看起来应该是这样的：

Original : This is my top secret message! Encoded : tspsm esos shit r imc g!eaTeye Decoded : This is my top secret message! Decoded correctly (1 true, 0 false) : 1

Encrypt

clear all;

Original='This is my top secret message!' a=randperm(length(Original)) Encoded=Original(a) p=[1:1:length(a)] x=[a;p]'

y=sortrows(x,1) Decorded=y(:,2)

s=Encoded(Decorded) disp(['Original:',Original]); disp(['Encoded:',Encoded]); disp(['Decoded:',s]);

correct=strcmp(Original,s);

disp(['Decoded correctly (1 true, 0 false)：',num2str(correct)])

第二题

目的：

3） 函数的使用以及数据可视化。

4） 题目要求要比题目一灵活，根据自己所想来完成可以实现题目要求的函数。

题目说明：所有的题目代码写入相应的脚本中，如果题目没有要求特定的脚本文件名，自行定义一个即可。

1.semilog plot：

在过去的5年中，某一个班级的学生数量分别为15，25，55，115，144。看起来，学生数量随时间是呈指数分布的。为了验证这个想法，使用semilogy对5年的学生数量进行绘制，使用xlabel，ylabel，title对图进行标注，曲线为虚线，红色，数据点使用方块为marker，marker颜色为绿色。请绘制该图，我想如果想要看到所有5个点的marker还需要使用xlim对x轴做标度。如果果真呈指数分布，那么结果通过semilogy的方式绘制出来看起来应该是呈线性的。

参考结果图：

103growth trendnumber of students1021010123years456

clear all; x=1:1:5;

y=[15,25,55,115,144];

semilogy(x,y,'rs--','MarkerFaceColor','g'); xlabel('years');

ylabel('number of students'); title('growth trend'); xlim([0 6]); 2.bar graph：

生成5个随机数的向量，使用bar函数绘制柱状图，要求bar为红色。结果图参考如下：

Bar graph of 5 random values10.90.80.70.60.50.40.30.20.1012345

clear all; x=1:1:5; y=rand(1,5); bar(x,y,'r');

title('Bar graph of 5 random values'); 3.插入和面绘制：

A．使用rand做一个5×5的随机矩阵z0；

B．使用meshgrid以及1：5的向量生成x0和y0，meshgrid的两个输入都使用1：5的向量； C．使用meshgrid以及1：0.1：5的向量生成x1和y1，两个输入都使用1：0.1：5的向量 D．使用interp2在x1和y1位置对x0，y0，z0进行插值得到z1； E．使用surf绘制z1，colormap设置为hsv，shading设置为interp； F．使用hold on保持figure，使用contour绘制等高线图； G．添加colorbar 结果参考下图：

10.910.80.60.40.20604020 004020600.80.70.60.50.40.30.20.10

clear all; z0=rand(5,5); x=1:5;

[x0,y0]=meshgrid(x,x); y=1:0.1:5;

[x1,y1]=meshgrid(y,y); z1=interp2(x0,y0,z0,x1,y1); surf(y,y,z1); colormap(hsv); shading(interp); hold on;

contour(y,y,z1); colorbar; hold off;

4.使用find：

写一个函数，函数声明如下：ind=findNearest(x，desiredVal)。x是一个向量，deisredVal为一个标量，该函数实现的功能为在x中找到与desiredVal最接近的数所在的位置：ind。如果x中有多个值接近desiredVal，则返回多个位置。使用的函数参考：abs，min和find。

函数运行结果可参考下例：

clear all; x=1:5;

desiredVal=2.5;

a=abs(x(:)-desiredVal) ind=min(find(a(:)<1)) 5.循环和流程控制

写一个函数叫loopTest(N)，该函数检测1：N的每一个数是否被2，3整除，如果被2整除则输出‘n is divisible by 2’，被3整除则输出’n is divisible by 3’，也可能都不能整除也能都不能整除则输出’n is NOT divisible by 2 AND 3.’ ‘n is divisible by 2 AND 3.’参考函数为for，mod，if，else，elseif。运行结果可参考下图：

function loopTest(6) for n=6

if mod(n,2)==0&&mod(n,3)~=0

disp(['',num2str(n) ,'n is divisible by 2'])

else if mod(n,3)==0&&mod(n,2)~=0 disp(['',num2str(n),'n is divisible by 2'])

else if mod(n,3)~=0&&mod(n,2)~=0 disp(['',num2str(n) ,'n is NOT divisible by 2 AND 3']) else

disp(['',num2str(n),'n is divisible by 2 AND 3']) end end end end

6.平滑器：

声明函数smoothed=rectFilt(x,width)，该函数将数据x以特定的窗宽范围内进行平均滤波。例如如果width=5，那么x中第n个点的值的平滑算法为mean（x(n-2:n+2)）。但是要注意的是在边缘会有泄漏的问题，上例中也就是在n<3和n>length(x)-2时会出问题。我想这几个边缘的点我们可以考虑保持原数据不动。

A． 注意经过平滑后得到的数据与原数据点数应相同； B． 为左右平均时对称，我们在输入参数width时应选择奇数； C． 确保边缘的情况不会出问题；

D． 提示使用循环和mean函数进行平滑运算，如果你对卷积很熟悉那么也可以使用conv

函数进行卷积运算；

E． 函数写完后，载入noisyData.mat数据进行测试，如果载入数据失败，那么使用xlsread

函数读取noisyData.xls数据也是一样的。 我使用width为11，参考图如下：

Smoothing Demonstration1.5Original DataSmoothed1 0.5Data value0-0.5-1-1.5 01020304050Index60708090100

x=xlsread('noisyData.xls'); width=11;

smoothed=rectFilt(x,width); figure; plot(x,'--'); hold on;

plot(smoothed,'r'); hold off;

xlabel('Index');

ylabel('Data value');

title('Smoothing Demonstration'); legend('Original Data','Smoothed'); 7.绘制圆的函数

我想绘制一个圆，matlab没有这个函数，没关系我们可以自己来写。编辑一个函数，[x,y]=getCircle(center,r)，其中x，y返回的是一个圆的x，y一一对应的值，圆心为center，center应该是一个2个元素的向量，第一个值为x坐标，第二个为y坐标，r为圆半径。最终我们利用返回的x和y只需使用plot函数就可以绘制一个圆。 提示：假如圆心为（0，0），那么半径为1的圆可以按照以下方程求解x，y坐标x(t)=cos(t), y(t)=sin(t)。其中t为[-2*pi:2*pi]。因此，我们可以在以上坐标基础上做任意的尺度变化（半径）和坐标移动（圆心坐标）。 完成以上函数后，我们可以用这个函数来做一些有意思的事情，请利用上面写出的函数完成下面两个练习：

A． 编写一个脚本，名为concentric.m，利用前面的圆函数，结合plot和hold on命令绘制5

个同心圆，结果参考下图：

1050-5-10-15-10-5051015

B． 编写一个脚本，名为olympic.m，利用上面的圆函数，绘制olympic的logo，结果如下

图：

0.80.60.40.20-0.2-0.4-0.6-0.8-1-1.2-1-0.500.511.5

程序(1)function [x,y]=getCircle(center,r) t=-2*pi:0.01:2*pi for i=1:length(t) a=t(i);

x(i)=r*cos(a)+center(1); y(i)=r*sin(a)+center(2); end

程序(2)r=0.5; center=[1,0];

[x,y]=getCircle(center,r) plot(x,y,'r','linewidth',3) hold on; r=0.5;

center=[-0.5,-0.5]; [x,y]=getCircle(center,r) plot(x,y,'y','linewidth',3) hold on; r=0.5;

center=[0.5,-0.5];

[x,y]=getCircle(center,r) plot(x,y,'g','linewidth',3) hold on; r=0.5;

center=[0,0];

[x,y]=getCircle(center,r) plot(x,y,'k','linewidth',3) hold on; r=0.5;

center=[-1,0];

[x,y]=getCircle(center,r) plot(x,y,'linewidth',3) xlim([-1.5 1.5]); ylim([-1.4 1]); hold off;

第三题

目的：

5） 通过自己写函数解决一些常见问题。

题目说明：所有的题目代码写入相应的脚本中，如果题目没有要求特定的脚本文件名，自行定义一个即可。

1.线性方程求解：

使用“\\”，求解下面的线性系统，计算并显示误差向量。

3a?6b?4c?1a?5b?27b?7c?3

clear all;

A=[3,6,4;1,5,0;0,7,7]; B=[1;2;3]; x=A\\B a=x(1) b=x(2)

c=x(3) 2.数值积分：

5使用trapz或者quad计算积分xe?x/3dx，我们知道根据分部积分法可以计算这个积分

?0得到?3xe?x/350?5/3那么请比较实用trapz计算得到的结果和直接使用|?9e?x/3|5??24e?9，0积分后数值计算的结果差异，并显示在命令行窗口中。

clear all; x=0:0.01:5; y=x.*exp(-x/3); z=trapz(x,y)

%%%%%%%%%%%%%%%%% syms a

int(a.*exp(-a/3),0,5)

3.矩阵求逆：

计算矩阵??12??的逆，使用inv函数。然后将结果乘以原始矩阵，检查结果是否为单34??位阵。

clear all; a=[1 2;3 4] b=a*inv(a)

4.多项式拟合：

使用xlsread函数分别载入数据randomDatax.xls和randomDatay.xls，使用ployfit函数对该两对变量分别做1，2，3，4，5阶拟合，为了得到比较好的拟合效果，使用polyfit的中心与标尺模式，具体使用方法为[p,S,mu] = polyfit(x,y,n)，它会返回三个参数，

接下来使用对应的[y,delta] = polyval(p,x,S,mu)来根据x对y进行拟合。最后绘制为图，原始数据为黑色圆点，其他5个拟合为不同颜色，对x，y轴做label，对整个图做title，另外加入legend以注明不同曲线的含义。

运行结果可参考下例：

Plynomial fits to random data0.20.150.10.050-0.05-0.1-0.15-0.2 100OriginalOrder 1Order 2Order 3Order 4Order 5110120130140150x160170180190200 y

clear all;

x=xlsread('randomDatax.xls'); y=xlsread('randomDatay.xls'); plot(x,y,'k.') hold on;

[p S mu]=polyfit(x,y,1);

[m delta]=polyval(p,x,S,mu); plot(x,m,'linewidth',2); hold on;

[p S mu]=polyfit(x,y,2);

[m delta]=polyval(p,x,S,mu); plot(x,m,'g','linewidth',2); hold on;

[p S mu]=polyfit(x,y,3);

[m delta]=polyval(p,x,S,mu); plot(x,m,'r','linewidth',2); hold on;

[p S mu]=polyfit(x,y,4);

[m delta]=polyval(p,x,S,mu); plot(x,m,'y','linewidth',2); hold on;

[p S mu]=polyfit(x,y,5);

[m delta]=polyval(p,x,S,mu); plot(x,m,'m','linewidth',2); hold off; xlabel('x'); ylabel('y');

title('Plynomialfits to randomdata');

legend('Original','Order 1','Order 2','Order 3','Order 4','Order 5',4); 5.神经脉冲传导过程Hodgkin-Huxley方程仿真

你需要先写一个ODE方程去描述神经脉冲。方程的描述方式我们参考了Hodgkin和Huxley1952年提出的模型。他们也因此获得了Noble奖。这个模型的思路是，神经膜有对电压敏感的闸门，随着膜电压的改变，这个闸门会打开或者关闭以形成离子通道。一旦闸门打开，离子就可以在他们之间流动从而影响膜电压。这个模型是非线性的，我们只能使用数值的方法求解。这个题目稍复杂，耐心按照下面一步步进行：

A． 首先下载HH.zip压缩包，解压后发现里面有6个m文件，分别是alphah.m，alpham.m，

alphan.m，betan.m，betam.m和betah.m。这些方程用来描述了闸门开转台alpha（V）和闭状态beta（V）与电压的关系。h,m,n是三个闸门。 B． 接下来按照下面的式子写ODE方程。

dn/dt?(1?n)?n(V)?n?n(V)dm/dt?(1?m)?m(V)?m?m(V)dh/dt?(1?h)?h(V)?h?h(V)dV/dt??1(Gkn4(V?Ek)?GNam3h(V?ENa)?GL(V?EL))C

其中C是膜电容，G是电导系数，E是K，Na，L通道的逆转电位。

C?1Gk?36GNa?120令，GL?0.3

Ek??72ENa?55EL??49.4有了这些常量以及上述微分方程，我们就可以写ODE方程了。 C． 编写一个脚本叫HH.m，我们使用ode45函数来求解微分方程，并且画出膜电压曲线。

首先我们定义时间尺度为0-20ms，我们的方程单位其实就是ms，初始状态为n=0.5;m=0.5,h=0.5;V=-60。计算结果应该为一个21×4的矩阵，其中每一列就是对应我们解出的n（t），m（t），h（t），和V（t）。我们这里绘制V（t）曲线，查看电压变化变至平稳的状态。参考图如下：

Approaching Steady State4020Transmembrane Voltage (mV)0-20-40-60-80024681012Time(ms)14161820

D． 更进一步，我们来考察一个神经元放电的情况，神经元只有在膜电压超过一定的电压阈

值才被认为有活动，为了找到这些阈值，我们求解这个模型10次，每次都用计算得到的结果作为下一次的初始解，不同的是，V在每次求解时都在上一次基础上加1，2，3，4，。。。10.求解完成后，我们找到每次求解结果的最大值，如果最大值大于0，则绘制为红色曲线，如果小于0则绘制为黑色。那么我们就可以从曲线的颜色上判断一个神经元是否发放了。参考图如下：

Threshold Behavior6040Transmembrane Voltage (mV)200-20-40-60-80024681012Time(ms)14161820

function dy=HHODEfun(t,y) %n,m,h,v C=1; Gk=36; Gna=120; Gl=0.3; Ek=-72; Ena=55; El=-49.4;

an=alphan(y(4)); bn=betan(y(4)); am=alpham(y(4)); bm=betam(y(4)); ah=alphah(y(4)); bh=betah(y(4)); dy=zeros(4,1);

dy(1)=(1-y(1))*an-y(1)*bn; dy(2)=(1-y(2))*am-y(2)*bm; dy(3)=(1-y(3))*ah-y(3)*bh;

dy(4)=-1/C*(Gk*y(1)^4*(y(4)-Ek)+Gna*y(2)^3*y(3)*(y(4)-Ena)+Gl*(y(4)-El));

figure(1)

[t,x]=ode45(@HHODEfun,0:0.001:20,[0.5 0.5 0.5 -60]);

plot(t,x(:,4))

xlabel('Time(ms)');

ylabel('Transmembrane Voltage (mV)'); title('Approaching Steady State');

figure(2) for i=1:10

[t,x]=ode45(@HHODEfun,0:0.001:20,[x(end,1) x(end,2) x(end,3) x(end,4)+i]); if max(x(:,4))>0

plot(t,x(:,4),'r','lineWidth',2); else

plot(t,x(:,4),'k','lineWidth',2); end hold on; end

xlabel('Time(ms)');

ylabel('Transmembrane Voltage (mV)'); title('Threshold Behavior');

6.Julia分形集：

这个练习我们试着产生Julia sets，关于Julia sets，对Julia sets感兴趣的同学可以参看http://en.wikipedia.org/wiki/Julia_Sets。其定义，给定两个复数c和z0，那么定义如下的递归：

2zn?zn?1?c，这个动态系统我们成为二次映射。给定c和z0的话，上述的递归方程会

给出一系列的z（k）值。不同的c和z0产生的结果当然是不一样的。这实际上也可以看作是一个混沌系统。我们下来来做这个系统，并且采用imagesc绘制看看这个美丽的分形图案。 A． 声明一个函数，julia（zMax，c，k，v）； B． 其中c是上述递归方程中的c，k给出了迭代次数，v指定了点数的规模，zMax为

了定义值的范围； C． 为了确定递归不会陷入无限循环，可以证明当z的模不大于2时，可以避免死循环，

因此在函数中应该声明一个常量为2；

D． 使用linspace，ones生成变量A，A的尺度为v×v的矩阵，其中的值在-zMax和

zMax之间。注意A是一个复数矩阵，复数的实部和虚部都应在-zMax和zMax之间；

E． 开始迭代，次数为k，在迭代中，每次根据上述递归方程得到新的A，如果A的模

不大于2，则保留，如此反复循环，将所有得到的新的模不大于2的A叠加在一起，最后是叠加完毕的v×v的矩阵。

F． 使用imagesc（0.1×atan（叠加后的矩阵））来显示你的结果，为了好看，可以另外

输入axis xy。

参考图如下：julia(1,-0.297491+i*0.641051,100,500)

5004504003503002502001501005050100150200250300350400450500

另，julia(.35,-.297491+i*0.641051,100,250)，c保持不变，点数变少，值范围变小，如下图：

2502001501005050100150200250

M=julia(1,-0.297491+i*0.641051,100,500); figure(1);

x=linspace(0,500,500); imagesc(x,x,atan(0,1*M)); axis xy

M=julia(.35,-.297491+i*0.641051,100,250); figure(2);

x=linspace(0,250,250); imagesc(x,x,atan(0,1*M)); axis xy

function n=escapeVelocity(z0,c,k,v) n=0; z=z0;

while abs(z)<=2&n

第四题

目的：

6） 练习高级功能，如数据结构，图像，动画等。

题目说明：所有的题目代码写入相应的脚本中，如果题目没有要求特定的脚本文件名，自行定义一个即可。

1.随机变量生成：

生成一个500个随机数的行向量（使用randn函数），通过计算使之变成一个均值为2，方差为5的随机数向量。之后，使用mean和std函数验证你所生成的随机数向量是否均值和方差接近2和5。

% normalRand

% homework 4 problem 1, making a vector of 500 random numbers that are % normally distributed with mean 2 and standard deviation 5 r=randn(500,1)*5+2;

% verify that the mean and std are close to 2 and 5 disp(['Mean: ' num2str(mean(r))]);

disp(['Standard Deviation: ' num2str(std(r))]);

2.硬币投掷试验：

建立一个脚本命名为coinTest.m，我们来模拟投掷硬币的正面概率问题，我们知道，如果硬币质量均匀，投掷次数足够多，那么两面出现的概率应该为50%，因此首先我们应该设定足够多的投掷次数，例如5000次。记录投出正面的次数，计算正面出现的概率，然后绘制曲线。每次投掷的结果可以用rand和round来模拟，假设1为正面，0则为背面。每次

投掷的结果累积和概率计算当然可以使用一个循环来完成。更好的办法是使用函数cumsum。另外把期望的0.5那条线也画出来以作对比。完成之后多运行几次，看看随着投掷次数的增加是不是正面出现的概率接近0.5。你也可以设置更多的投掷次数，例如下图是我做了10000次的结果。

% coinTest

% generate 1000 random numbers r=rand(5000,1);

% round them so 1 means heads r=round(r);

% count up the number of heads and divide by the count count=1:length(r);

headProb=cumsum(r)./count'; % plot the probability of heads figure

plot(count,headProb,'k','linewidth',1.5); hold on

plot([1 length(r)],[.5 .5],'r--');

title('Sample Probability of Heads in n flips of a simulated coin'); xlabel('Number of coin flips'); ylabel('Probability of heads');

legend('Sample Probability','Fair coin');

sample probability of Heads in n flips of a simulated coin1sample probabilityFair coin 0.90.8probability of heads0.70.60.50.4 0100020003000400050006000Number of coin flips70008000900010000

3.cell练习：

通常cell被用来存储字符串，因为对于每个位置的字符串可以使用不同的长度。 1）建立一个3×3的cell，第一列赋值为‘Joe’，‘Sarah’，‘Pat’。第二列包含了他们的last name：‘Smith’，‘Brown’，‘Jackson’。第三列则是他们各自的收入状况：‘$30,000’, ’$150,000’, ‘$120,000’。赋值后显示结果在命令行查看是否正确。

2）Sarah结婚了将她的last name换成了‘Meyers‘，对之前生成的cell进行修改并显示查看。

3）Pat涨了$50，000的工资，修改她的收入，并显示查看。 结果参考图如下： cellPro =

'Joe' 'Smith' '$30,000' 'Sarah' 'Brown' '$150,000' 'Pat' 'Jackson' '$120,000'

cellPro =

'Joe' 'Smith' '$30,000' 'Sarah' 'Meyers' '$150,000' 'Pat' 'Jackson' '$120,000'

cellPro =

'Joe' 'Smith' '$30,000' 'Sarah' 'Meyers' '$150,000' 'Pat' 'Jackson' '$170,000'

% cellProblem

% initializes a cell, puts some values in it, and changes some of them % part a

% initialize the cell c=cell(3,3); % add elements c{1,1}='Joe'; c{2,1}='Sarah'; c{3,1}='Pat'; c{1,2}='Smith'; c{2,2}='Brown'; c{3,2}='Jackson'; c{1,3}=30000; c{2,3}=150000; c{3,3}=120000; % display the cell disp(c)

% part b

% change sarah's name to Meyers c{2,2}='Meyers'; disp(c) % part c

% add 50000 to pat's salary c{3,3}=c{3,3}+50000; disp(c)

4.使用structs：

Structs在处理多样数据时非常有用。例如在命令行输入a=dir，可以看到返回的a就是一个structs，其中包括了很多域，例如name，bytes，isdir等。

1）使用a.name, a.bytes等查看当前路径下文件的情况；

2）编写一个循环代码，遍历a中每一个文件，如果这个文件不是一个文件夹，则输出如下文本“File name contains X bytes”。

运行结果可参考下例：

HW4.doc contains 127488 bytes Pro1.m contains 58 bytes cellPrc.asv contains 161 bytes cellPrc.m contains 198 bytes coinTest.asv contains 374 bytes coinTest.m contains 382 bytes structPrc.asv contains 79 bytes structPrc.m contains 151 bytes ~$HW4.doc contains 162 bytes

~WRL0001.tmp contains 128000 bytes

% displayDir %

% displays a list of the files in the current directory and the size of % each file in bytes function displayDir

% get current directory info a=dir;

% loop through it and display each file name and size for n=1:length(a)

% if it's not a directory, display the information if ~a(n).isdir

disp(['File ' a(n).name ' contains ' num2str(a(n).bytes) ' bytes.']); end end

5.Handels使用：

使用handels来完成一个如上图的图样来： 1）建立一个脚本命名为handelsPractice.m；

2）首先生成一个变量x，0：0.01：2pi，生成y=sin（x）； 3）使用figure打开一个绘图界面，使用plot（x，y，‘r’），绘制正弦曲线‘ 4）使用xlim设置x轴范围为0～2pi；

5）使用set和gca（get current axis）来设置xtick属性，值设为[0 pi 2*pi],设置xticklabel属性，值为{‘0’,’1’,’2’};

6）使用set和gca设置ytick属性，值为-1：0.5：1； 7）使用grid on打开网格；

8）使用set和gca设置ycolor属性为green，xcolor属性为cyan，另外设置color属性值为black；

9）使用set和gcf（注意这里是gcf，get current figure），设置figure的color属性为[0.3 0.3 0.3]；

10）使用title添加标题“One sine wave from 0 to 2π”，fontsize属性为14，fontweight属性为bold，color属性为white，为了显示为π，使用\\pi来正确显示这些数学符号；

11）设置xlabel和ylabel，fontsize属性为12，color属性分别为cyan和green； 12）如果你想粘贴到word中保持我们调成的这个样式，在copy figure之前，在copy options中figure background color中选择use figure color。至此完成。

% handlesPractice % plot the function x=linspace(0,2*pi,100); y=sin(x); figure;

plot(x,y,'r'); xlim([0 2*pi]); % set the ticks

set(gca,'xtick',[0 pi 2*pi],'xticklabel',{'0','1','2'}); set(gca,'ytick',[-1:.5:1]);

% turn on grid and set axis and background colors grid on;

set(gca,'ycolor','g') set(gca,'xcolor','c'); set(gca,'color','k');

set(gcf,'color',[.3 .3 .3]); % add the labels

title('One sine wave from 0 to

2\\pi','fontsize',14,'fontweight','bold','color','w')

xlabel('X values (in terms of \\pi)','fontsize',12,'color','c'); ylabel('Sin(X)','fontsize',12,'color','g');

6.图片读取与绘制：

使用imread读入12110015.jpg（我家两只猫的照片）。然后绘制原图与其R，G，B分量的图，共四张图。当你用imread读入图片后会发现这是一个三维矩阵，每一页的维度表示像素大小，三页分别表示R，G，B的颜色编码。如果你只要绘制红色，可以将后两页的值置为0来进行绘制。使用subplot，绘制为2×2的样式。结果参考下图：

500100015002000500100015002000100020003000100020003000500100015002000500100015002000100020003000100020003000 % im=displayRGB(filename)

% reads in the jpg file in filename (filename should include the % extension), displays the original image as well as each color layer in a

% figure separately, and returns this composite image in im. the original % image is resampled so the largest final dimension is 800 pixels and the % aspect ratio is preserved. x=imread('12110015.jpg'); subplot(2,2,1),imshow(x); a=x(:,:,1); aa=x(:,:,2); aaa=x(:,:,3);

r=cat(3,a,zeros(size(a)),zeros(size(a))); subplot(2,2,2),imagesc(r);

g=cat(3,zeros(size(aa)),aa,zeros(size(aa))); subplot(2,2,3),imagesc(g);

b=cat(3,zeros(size(aaa)),zeros(size(aaa)),aaa); subplot(2,2,4),imagesc(b); end

7.动画模拟Brown运动：

声明一个函数叫brown2D（N）。这个函数只有一个输入参数为N，表示用来模拟Brown运动的点数。每一个点初始位置为（0，0），plot所有的点使用‘.‘作为marker。设置xlim和ylim为-1～1，axis设置为square。做1000次模拟，即写一个1000次的循环。每次更新

点的x，y坐标值，更新的位置使用randn生成，方差为0.005。加入N为100，那么对于每一个x，y都需要使用randn来更新其坐标值。更新的值等于此次生成的随机变量值加上上一次的坐标值，在播放动画时，使用pause来停顿0.01秒。最好使用set来设置xdata和ydata属性的方式来更新x，y的坐标值。运行程序会看到点从中心开始向四周方向随机漂移。结果可参考我制作的avi视频，brown2D.avi。

% brown2D(N) %

% simulates the Brownian motion of N points. all N points begin in the same

% location and diffuse away from it as time passes. The animation runs for

% about 10 seconds function brown2D(N) % make a figure figure; axis square

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

% make the x and y positions of the points x=zeros(N,1); y=zeros(N,1); h=plot(x,y,'.k');

% run a loop to display the points for n=1:1000

% update the positions x=x+randn(size(x))*.005; y=y+randn(size(y))*.005; % plot them

set(h,'xdata',x,'ydata',y); axis([-1 1 -1 1]) pause(0.01); end

8.Julia分形集动画制作：

第三次做的Julia函数可以用来绘制julia分形图，现在改成一个动画，我们只需要在迭代过程中，对每次更新的A矩阵都做绘制，使用drawnow，这样会把每一次更新的A都进行更新绘制就形成了一个很酷的动画。这里不使用pause是因为julia分形的迭代计算耗费较长的时间，因此我们无法定义pause多长时间，而使用drawnow则可以无需等待，很快的更新结果。程序运行结果参考我制作的avi视频，JuliaAni.avi。

% juliaAnimation(zMax,c,N) %

% displays an animation of a Julia Set %

% inputs:

% zMax - the complex numbers for which escape velocities are computed will

% be 500 values between -zMax and zMax for both real and imaginary % part

% c - the complex number c, which is a parameter in the Julia set % N - the maximum number of iterations to use when computing escape % velocities; also equals the number of frames in the animation %

% EXAMPLES

% juliaAnimation(1,-.4+i*.6,75); % juliaAnimation(1,-.8+i*.156,100);

% juliaAnimation(.75,-.297491+i*.641051,150); function juliaAnimation(zMax,c,N)

% make the complex grid, with real and imaginary part going from -zMax to % zMax

res=500; % resolution, number of points on a side of the grid temp=linspace(-zMax,zMax,res); [R,I]=meshgrid(temp,temp); Z=R+1i*I;

% initialize the map to have 'escape velocities' all equal to N M=N*ones(size(Z));

% loop N times, each time generating the new Z and visualizing it figure for n=1:N

% calculate Zn for the current n Z=Z.^2+c;

% find where the modulus is >2 inds=find(abs(Z)>2);

% set these inds in M to contain the loop iteration M(inds)=n;

% set these inds in Z to be nan Z(inds)=nan;

imagesc(temp,temp,atan(0.1*M)); axis xy drawnow

end

本文来源：https://www.bwwdw.com/article/wyg7.html