matlab作业3

信息科学与工程学院
Matlab作业
**: ***
专业: 电子信息工程2班
学号: ********** 日期 2015 年 6 月 28 日
1.线性方程求解：
使用“\”，求解下面的线性系统，计算并显示误差向量。

3641
52
773
a b c a b b c ++=+=+= sy3_1.m
A=[3,6,4;1,5,0;0,7,7]; B=[1;2;3]; X=A\B w=A*X-B 运行结果：X = -0.5824 0.5165 -0.0879 w =
1.0e-15 *
-0.1110 0 0
2.数值积分：
使用trapz 或者quad 计算积分5
/3
x xe
dx -⎰
，我们知道根据分部积分法可以计算这个积分
得到/35/35
5/3003|9|249x x xe e
e -----=-+，那么请比较实用trapz 计算得到的结果和直接使用积分后数值计算的结果差异，并显示在命令行窗口中。

Sy3_2.m
x=0:0.0001:5; y=x.*exp(-x/3); z0=trapz(x,y) z1=-24*exp(-5/3)+9 w=z0-z1
运行结果：
z0 = 4.4670
z1 = 4.4670 w =
-9.3826e-10 3.矩阵求逆：
计算矩阵
12
34
⎡⎤
⎢⎥
⎣⎦
的逆，使用inv函数。

然后将结果乘以原始矩阵，检查结果是否为单
位阵。

Sy3_3.m
a=[1,2;3,4];
b=inv(a)
c=b*a
运行结果：
b =
-2.0000 1.0000
1.5000 -0.5000
c =
1.0000 0
0.0000 1.0000
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以注明不同曲线的含义。

Sy3_4.m
x=xlsread('randomDatax.xls');
y=xlsread('randomDatay.xls');
plot(x,y,'k.')
hold on
color='bgrym';
for n=1:5
[p,S,mu] = polyfit(x,y,n);
[y1,delta] = polyval(p,x,S,mu);
plot(x,y1,color(n),'linewidth',2)
hold on
end
xlabel('x');
ylabel('y');
title('plynomial fits to random data');
legend('Original','Order1','Order2','Order3','Order4','Order5');
运行结果：
100
110120130140
150160170180190200
-0.2-0.15-0.1-0.0500.050.10.150.2x
y
plynomial fits to random data
Original Order1Order2Order3Order4Order5
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 方程。

43/(1)()()/(1)()()/(1)()()1
/(()()())n n m m h h k k Na Na L L dn dt n V n V dm dt m V m V dh dt h V h V dV dt G n V E G m h V E G V E C
αβαβαβ=--=--=--=-
-+-+-
其中C 是膜电容，G 是电导系数，E 是K ，Na ，L 通道的逆转电位。

令，
1
36
120
0.3
72
55
49.4
k
Na
L
k
Na
L
C
G
G
G
E
E
E
=
=
=
=
=-
=
=-
有了这些常量以及上述微分方程，我们就可以写ODE方程了。

ODE.m
function dy = ODE(t,y)
C=1;
G_k=36;
G_Na=120;
G_L=0.3;
E_k=-72;
E_Na=55;
E_L=-49.4;
n=y(1);
m=y(2);
h=y(3);
V=y(4);
dy=zeros(4,1);
dy(1)=(1-n)*alphan(V)-n*betan(V);
dy(2)=(1-m)*alpham(V)-m*betam(V);
dy(3)=(1-h)*alphah(V)-h*betah(V);
dy(4)=-(1/C)*((G_k*n^4*(V-E_k))+G_Na*m^3*h*(V-E_Na)+G_L*(V-E_L)); end
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）曲线，查看电压变化变至平稳的状态。

D．更进一步，我们来考察一个神经元放电的情况，神经元只有在膜电压超过一定的电压阈值才被认为有活动，为了找到这些阈值，我们求解这个模型10次，每次都用计算得到的结果作为下一次的初始解，不同的是，V在每次求解时都在上一次基础上加1，2，3，4，。

10.求解完成后，我们找到每次求解结果的最大值，如果最大值大于0，则绘制为红色曲线，如果小于0则绘制为黑色。

那么我们就可以从曲线的颜色上判断一个神经元是否发放了。

：
HH.m
figure(1);
[t,y]=ode45(@ODE,[0 20],[0.5 0.5 0.5 -60]);
n=y(:,1);
m=y(:,2); h=y(:,3); V=y(:,4); plot(t,V);
xlabel('Time(ms)');
ylabel('Transmembrane(mV)');
title('Approaching Steady State'); figure(2); for i=1:10
[t,y]=ode45(@ODE,[0 20],[n(end) m(end) h(end) V(end)+i]); n=y(:,1); m=y(:,2); h=y(:,3); V=y(:,4); if (max(V)>0)
plot(t,V,'r','linewidth',2); hold on ;
else plot(t,V,'k','linewidth',2); hold on ; end end
xlabel('Time(ms)');
ylabel('Transmembrane(mV)'); title('Approaching Steady State'); hold off ; 运行结果：
2
4
6
8
101214
16
18
20
-80-60-40
-20
20
40
Time(ms)
T r a n s m e m b r a n e (m V )
Approaching Steady State
02468
101214161820
-80
-60
-40
-20
2040
60
Time(ms)
T r a n s m e m b r a n e (m V )
Approaching Steady State
6.Julia 分形集：
这个练习我们试着产生Julia sets ，关于Julia sets ，对Julia sets 感兴趣的同学可以参看/wiki/Julia_Sets 。

其定义，给定两个复数c 和0z ，那么定义如下的递归：
2
1n n z z c -=+，这个动态系统我们成为二次映射。

给定c 和0z 的话，上述的递归方程会
给出一系列的z （k ）值。

不同的c 和0z 产生的结果当然是不一样的。

这实际上也可以看作是一个混沌系统。

我们下来来做这个系统，并且采用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.m
function julia(Zmax,c,k,v)
d=linspace(-Zmax,Zmax,v);
A=ones(v,1)*d+i*(ones(v,1)*d)';
B=zeros(v,v);
for s=1:k
B=B+(abs(A)<=2);
A=A.*A+ones(v,v).*c;
end
imagesc(0.1*atan(B));
colormap(jet);
axis xy;
end
julia(1,-0.297491+i*0.641051,100,500)
运行结果：
500
450
400
350
300
250
200
150
100
50
50100150200250300350400450500
julia(.35,-.297491+i*0.641051,100,250)，c保持不变，点数变少，值范围变小运行结果：
250
200
150
100
50
50100150200250。