HH Model PDE (errors)

조회 수: 2 (최근 30일)
Hyeok Jun Lee
Hyeok Jun Lee 2019년 6월 27일
Hi I'm pretty new at MATLAB, and I'm trying to solve a system of PDEs involving the cable equation using pdepe function. However I get an error saying:
"Unable to meet integration tolerances without reducing the step size below the smallest value allowed" and that time integration failed. I was wondering if there is any way around this or if there is a way to estimate the PDEs recursively or through another method. Here is the code I have written for it:
function HH
m = 0;
x = linspace(-30,30,150);
t = linspace(0,200,50);
sol = pdepe(m,@pde_eqn,@pde_initial,@pde_bc,x,t);
V = sol(:,:,1);
%produce Figure 8.12B
figure(2)
set(gca,'fontsize',14)
hold on
plot(x, V(1,:),'k', 'linewidth', 2);
plot(x, V(50,:),'k--', 'linewidth', 2);
plot(x, V(150,:),'k:', 'linewidth', 2);
axis([0 15 -80 40])
xlabel('Position (cm)')
ylabel('Membrane voltage (mV)')
legend('t=0 msec', 't=50 msec', 't=150 msec')
str1(1) = {'B'};
text(-2.5,60,str1, 'Fontsize', 40)
end
function [c,b,s] = pde_eqn(x,t,u,DuDx)
E_Na = 52.4;
E_K = -72.1;
E_L = -49.2;
g_Na = 120;
g_K = 36;
g_L = 0.3;
c_m = 1.0;
Ri = 35;
Re = 20;
a_m = @(x)(0.1*(x+40)/(1-exp(-0.1*(x+40))));
b_m = @(x)(4*exp(-0.0556*(x+65)));
a_n = @(x)(0.01*(x+55)/(1-exp(-0.1*(x+55))));
b_n = @(x)(0.125*exp(-0.0125*(x+65)));
a_h = @(x)(0.07*exp(-0.05*(x+65)));
b_h = @(x)(1/(1+exp(-0.1*(x+35))));
c = [1;1;1;1];
b = [1/(c_m*(Ri+Re)); 0; 0; 0].*DuDx;
s = [(1/c_m)*(g_Na*u(2)^3*u(4)*(u(1)-E_Na)+g_K*u(3)^4*(u(1)-E_K)+g_L*(u(1)-E_L)); a_m(u(1))*(1-u(2))-b_m(u(1))*u(2); a_n(u(1))*(1-u(3))-b_n(u(1))*u(3); a_h(u(1))*(1-u(4))-b_h(u(1))*u(4)];
end
function value = pde_initial(x)
a_mi = @(y)(0.1*(y+40)/(1-exp(-0.1*(y+40))));
b_mi = @(y)(4*exp(-0.0556*(y+65)));
a_ni = @(y)(0.01*(y+55)/(1-exp(-0.1*(y+55))));
b_ni = @(y)(0.125*exp(-0.0125*(y+65)));
a_hi = @(y)(0.07*exp(-0.05*(y+65)));
b_hi = @(y)(1/(1+exp(-0.1*(y+35))));
value = [-65; a_mi(0)/(a_mi(0)+b_mi(0)); a_ni(0)/(a_ni(0)+b_ni(0)); a_hi(0)/(a_hi(0)+b_hi(0))];
end
function [pl, ql, pr, qr] = pde_bc(xl,ul,xr,ur,t)
pl=[0;0;0;0];
ql=[1;1;1;1];
pr=[0;0;0;0];
qr=[1;1;1;1];
end

답변 (0개)

카테고리

Help CenterFile Exchange에서 Boundary Conditions에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by