How to insert initial condition

Hello everybody, I trying to solve cable equation with PDE solver i would like to add punctually a temporal condition, but for beginning i'm trying to just set an array of values in the initial condition of my function. My equation looks like "(alpha)*dV/dt = (beta)*d²v/d²t - v" ;

if true
function Cable_transport_HH
m = 1 ; % cylinder
nx = 20 ; %spatial discretization
nt = 100 ; % temporal discretisation
x = linspace(0,500e-6,nx);%spatial space
t = linspace(0,20e-3,nt);% time space
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);% pde solver
u = sol(:,:,1);% solution
function [c,f,s] = pdex1pde(x,t,u,DuDx)
% constant
Cm = 1 ; % Membrane Capcitance mF/cm^2
Rm = 2500e-3 ; % Ohm.cm²
Ri = 70e-3 ;% Ohm.cm
d = 5e-4 ; % diameter cm 
alpha = pi*d*Cm ; 
beta = (pi*d^2)/(4*Ri) ;
gamma = (pi*d/Rm) ;
c = alpha;
f = beta*DuDx;
s = gamma*u;
function u0 = pdex1ic(x)
u0 = IT IS here that i should add my array of value for V(0,t)?? ;
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = 1; %%here i set arbitrary but my cable on the right is open and i haven't found how to set this 
ql = 1;
pr = 1;
qr = 1;
end

Is someone could help me to understand the PDE for solving this equation. Thanks in advance, Best regards, Antoine

[Torsten]

I don't see an x-variable in your equation - so why do you use "pdepe" instead of "bvp4c"?

Is "v" = "V" in your equation ?

Does d²v/d²t mean d²v/dt² ?

What are your boundary conditions ?

Best wishes

Torsten.

[antoine jury]

aha because i didn't knew this function i will try to look how to do with thanks

[antoine jury]

And sorry the equation was "alpha*dv/dt = beta*dv²/dx² -gamma*v", that's the equation with the x.

[Torsten]

I don't know your boundary conditions, but maybe already ODE45 suffices to solve your problem.

Best wishes

Torsten.

[antoine jury]

Thanks for your answer. My boundary conditions are simple i have some array of data which define the v(0,t) and the point v(L,t) is open.

[Torsten]

You need two boundary conditions to solve your problem.

If both of them are given either at t=0 or at t=L, you can use ODE45.

If one condition is given at t=0 and the other at t=L, you will have to use bvp4c.

Best wishes

Torsten.

[antoine jury]

i had found one question similar using pdede but the boundary is function time dependent.

    function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
    pl = T1( t );
    pr = T2;
    ql = 0;
    qr = 0;
function T = T1( t )
    T = ...

In my problem the value of V(0,t) is calculated from a system of ode.

• tau*dv/dt = -v + n(w-v)/Rf
• Cf*dw/dt = (v-w)/Rf +qI
• dm/dt = phi*[a*(1-m) - B(m)*m]
• dh/dt = phi*[a(h)*(1-h) - B(h)*h]
• dn/dt = phi*[a(n)*(1-n) - B(n)*n]

By this way I solve at one point (v(0,t)) my problem but than the previous equation ("alpha*dv/dt = beta*dv²/dx² -gamma*v"), describe the conduction of this input signal. If you could help me again to understand it's could be kind from you .

Best regards,

Antoine.

[Torsten]

Since you have an ordinary differential equation

tau*dv/dt = -v + n(w-v)/Rf

to define the boundary value for v at x=0, you can't use PDEPE.

You will have to discretize the expression d^2v/dx^2 in space and solve the resulting system of ordinary differential equations for v in the grid points in x-direction together with the five ODEs stated for v,w,m,h,n using ODE15S.

Look up "method-of-lines" for more details.

Best wishes

Torsten.

[antoine jury]

I already try by this way but i didn't succed to figure out how to solve my problem, i send you my script, i don't understand how the final matrix is construct.

clc; clear all; tic ;
%%Constants set
Cm = 1 ; tf = 20 ; I = 0.1 ; Cf = Cm ; 
ENa = 115 ; EK = -12 ; El = -3.8038 ; 
gbarNa = 180 ; gbarK = 18 ; gbarl = 0.05 ; 
mu = 5 ; q = 8 ;
Rm = 2500 ; Ri = 70 ; tau = Rm*Cm ;
l = 0.5e-4 ; df = 0.1e-4 ; d = 5e-4 ;
Rf = (4*Ri*l)/(pi*df^2) ; 
T = 37 ; Q10 = 3^((T-6.3)/10) ; k = (pi*d^2) / (4*Ri) ;
%%ODE45 Method
V = 0 ; % Initial Membrane voltage mV
W = 0 ; 
m = am(W) / ( am(W)+bm(W) ) ; % Initial m-value
n = an(W) / (an(W)+bn(W)) ; % Initial n-value
h = ah(W) / (ah(W)+bh(W)) ; % Initial h-value
y0 = [ V ; W ; n ; m ; h ] ; %Vector of initial conditions
tspan = [0,tf] ;
%Matlab's ode function
[time,V] = ode15s(@(t,y) BH_conduction(t,y), tspan, y0);
ODV = V(:,1) ; ODW = V(:,2) ; ODn = V(:,3) ; ODm = V(:,4) ; ODh = V(:,5) ;
clear V ;
%%Plots
%Plot the functions
figure
plot(time,ODV)
legend('ODE')
xlabel('Time (ms)')
ylabel('Voltage (mV)')
title('Voltage Change for Hodgkin-Huxley Model')
figure
plot(time,ODn,'b--',time,ODm,time,ODh,'r--');
ylabel('Gaining Variables')
xlabel('Time (ms)')
axis([0 tf 0 1])
legend('n','m','h');
toc;

And the function that I'm calling

function fval = BH_conduction(t,y)
Cm = 1 ; tf = 20 ; I = 0.1 ; Cf = Cm ; 
ENa = 115 ; EK = -12 ; El = -3.8038 ; 
gbarNa = 180 ; gbarK = 18 ; gbarl = 0.05 ; 
mu = 5 ; q = 8 ;
Rm = 2500 ; Ri = 70 ; tau = Rm*Cm ;
l = 0.5e-4 ; df = 0.1e-4 ; d = 5e-4 ;
Rf = (4*Ri*l)/(pi*df^2) ; 
T = 37 ; Q10 = 3^((T-6.3)/10) ; k = (pi*d^2) / (4*Ri) ;
% Values set to equal input values
V = y(1); W = y(2); n = y(3); m = y(4); h = y(5);
gNa=gbarNa*m^3*h; gK=gbarK*n^4; gl=gbarl;
INa=gNa*(V-ENa); IK=gK*(V-EK); Il=gl*(V-El);
% parameters 
alpha = 1/(pi*df*Cm) ; beta = (pi*df^2)/(4*Ri) ; gamma = (pi*df)/Rm ;
N = length(y)+1 ;
DyDt = zeros(length(y),N) ;
for i = 2:N
  DyDt = [  (alpha)*( beta*(V(1,i+1)-2*V(1,i)+V(1,i-1) ) - gamma*V(1,i) + (W-V(1,i))/Rf ) ; ...
            ((1/Cf)*(I-(INa+IK+Il))); ...
            an(W)*(1-n)-bn(W)*n; ...
            am(W)*(1-m)-bm(W)*m; ...
            ah(W)*(1-h)-bh(W)*h];
% extraction fval from Dv/Dt 
fval = DyDt(2:N) ;

Thanks in advance if you have time to check my script.