Error using pdepe function in matlab.

The code runs well with M as a numeric but when I change it to a range of values, it gives me a error. How can I rectify this?

Answer 1

[Torsten]

Call pdepe in a loop for M = M_array(1),M = M_array(2),...

Best wishes

Torsten.

[Brian]

Thanks Torsten. I have not understood what you mean. Could you please explain a bit more? Thanks again!

[Torsten]

function trial
M = [0,2.5950,4.5412,5.8387,6.4874,6.4874,5.8387,4.5412,2.5950,0];
% Define the solution mesh here
x = linspace(0,1,10); 
t = linspace(0,1,10); 
options = odeset('RelTol',2e+02);
%--------------------------------------------------------------------------
m = 0; 
for k=1:numel(M)
  M_actual = M(k);
  sol{k} = pdepe(m, @(x,t,u,DuDx)fpde(x,t,u,DuDx,M_actual), @fic, @fbc,x,t, options);
end
...

function [c,f,s] = fpde(x,t,u,DuDx,M) 
% %---------------------------------------------------------------------
ce = 4.0; 
cs = 125; 
lambdaphi = 7.91*10^(-18);
Dphi = 5*10^(-7); 
dk = 8.64*10^(4); 
ds = 0.25; 
de = 2.16*10^(4); 
lambdam = 2.0;
phic = 200; 
u(1) =10^(-3);
u(2) =9*u(1);
Ku = 0.8;
%--------------------------------------------------------------------------
u(3) = 1;
lambda0 = lambdaphi.*u(3).*(1+(u(1)/u(2)));
gamma = de/dk;
delta = ds/dk;
eta = (lambdaphi/lambda0)*phic;
m0= lambdam/Ku;
%--------------------------------------------------------------------------
lambdau = 4.0; 
Du = 5*10^-5;
Km = 0.28;     
Dm = 10^-2;
beta = Km/lambdau;
%--------------------------------------------------------------------------
Me = ce./m0;
Ms = cs./m0;
z1 = heaviside(M-Me);
z2 = heaviside(M-Ms);
z3 = heaviside(eta-u(3));
z11 =(-gamma .*u(1).*z1)- (z3.*u(1))';
z22 =(-gamma .*u(1).*z1)- (z3.*u(1))';
z33 = u(1) - (u(2)+u(1)).* u(3)';
% Left hand of the pdes.
%--------------------------------------------------------------------------
c = [1;1;1]; 
s = [z11; z22; z33];
f = [0;
     0;
     DuDx(3)];
end

...

[Brian]

It works. Thanks alot, Torsten!

[Brian]

Dear Torsen,

The code did run but did not show a change in the illustration. I would like to plot M(x) in ee and ss. M(x) is the myco.csv file and it has two columns. I expect M to iterate on two different values from each row. M(x) is illustrated as the last graph plotted. I expect to see the straight line in the middle part of the first two graphs distorted such as a bulge(Because that how M(x) looks like). I have attached the csv file and the code. I hope I am clear enough. Thanks a lot again. Cheers!

[Torsten]

You never told us that M is a function of x. As this seems to be the case, my suggested solution is wrong.

Insert xM in the following code where xM is the x-value where M(x) is attained. I guess these are the values from your Excel-sheet.

function trial
xM = [...];
M = [0,2.5950,4.5412,5.8387,6.4874,6.4874,5.8387,4.5412,2.5950,0];
% Define the solution mesh here
x = linspace(0,1,10); 
t = linspace(0,1,10); 
options = odeset('RelTol',2e+02);
%--------------------------------------------------------------------------
m = 0; 
sol = pdepe(m, @(x,t,u,DuDx)fpde(x,t,u,DuDx,xM,M), @fic, @fbc,x,t, options);

...

function [c,f,s] = fpde(x,t,u,DuDx,xM_array,M_array) 
% %---------------------------------------------------------------------
M=interp1(xM_array,M_array,x);
ce = 4.0; 
cs = 125; 
lambdaphi = 7.91*10^(-18);
Dphi = 5*10^(-7); 
dk = 8.64*10^(4); 
ds = 0.25; 
de = 2.16*10^(4); 
lambdam = 2.0;
phic = 200; 
u(1) =10^(-3);
u(2) =9*u(1);
Ku = 0.8;
%--------------------------------------------------------------------------
u(3) = 1;
lambda0 = lambdaphi.*u(3).*(1+(u(1)/u(2)));
gamma = de/dk;
delta = ds/dk;
eta = (lambdaphi/lambda0)*phic;
m0= lambdam/Ku;
%--------------------------------------------------------------------------
lambdau = 4.0; 
Du = 5*10^-5;
Km = 0.28;     
Dm = 10^-2;
beta = Km/lambdau;
%--------------------------------------------------------------------------
Me = ce./m0;
Ms = cs./m0;
z1 = heaviside(M-Me);
z2 = heaviside(M-Ms);
z3 = heaviside(eta-u(3));
z11 =(-gamma .*u(1).*z1)- (z3.*u(1))';
z22 =(-gamma .*u(1).*z1)- (z3.*u(1))';
z33 = u(1) - (u(2)+u(1)).* u(3)';
% Left hand of the pdes.
%--------------------------------------------------------------------------
c = [1;1;1]; 
s = [z11; z22; z33];
f = [0;
     0;
     DuDx(3)];
end

...

Best wishes

Torsten.

[Brian]

Perfect! Thanks!

[Brian]

I am sorry I am asking again and again! I made several adjustments in my code and got stuck again! I realized the Heaviside function couldn't do what I wanted it to do so I wrote functions for the Heaviside function and I ended up with s that is of a different dimension with f and c. Can someone help me figure that out. Thanks a bunch!

[Walter Roberson]

Your code does not explicitly define a value for eta == cee exactly. If the final value happens to be an exact match, then your hvf will come up one element short of the length it needs to be.

Consider

clear hvf
hvf([true false false]) = 3
hvf([false true false]) = 5

this does not give you a vector of length 3 on output with the third element being 0: this gives you a vector of length 2, because trailing false in a logical selector never extend the length of the variable on assignment.

You should be initializing

hvf = zeros(size(cee));

and that would ensure that any value not assigned to will be 0.

[Walter Roberson]

Torsten had suggested

z1 = heaviside(M-Me);
z2 = heaviside(M-Ms);
z3 = heaviside(eta-u(3));

where M was the interpolated value based upon x. You replaced the heaviside with calls to Untitl that reads the entire .cvs file and uses all of it instead of using the interpolated value -- and without even differentiating between the two columns of the file at that.

Do not read the entire file there. Pass in the current M value to your Untitl function and use that.

I notice, by the way, that your heav function calculates a constant eta and compares that to the input parameter, and that your weakly-named Untitl just does comparisons, acting as a heaviside function. I would suggest to you that it would make more sense to have a function that calculated eta once and return that, and then to have a heaviside function that did nothing but heaviside, and pass in the appropriate difference between eta and whatever. Or don't bother with an explicit heaviside and instead code

z1 = M > Me;
z2 = M > Ms;
z3 = eta() > u(3);

where

function eta_val = eta()
  phi =100;
  lambdaphi = 4.9*10^(-3);
  phic = 60;
  lambda0 = 10* lambdaphi.*phi;
  eta_val = (lambdaphi/lambda0)*phic;

[Brian]

Thanks Walter! This helps.