HH Model PDE (errors)
조회 수: 2 (최근 30일)
이전 댓글 표시
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
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Boundary Conditions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!