[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);