Error using pdepe function in matlab.

Question

Brian 2018-2-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/381755-error-using-pdepe-function-in-matlab

评论： Brian 2018-2-20

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2018-2-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/381755-error-using-pdepe-function-in-matlab#answer_304338

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

Best wishes

Torsten.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Torsten 2018-2-9

编辑：Torsten 2018-2-9

在 MATLAB Online 中打开

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
...

Torsten 2018-2-16

编辑：Torsten 2018-2-16

在 MATLAB Online 中打开

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.

Walter Roberson 2018-2-20

在 MATLAB Online 中打开

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 2018-2-20

Thanks Walter! This helps.

请先登录，再进行评论。

Error using pdepe function in matlab.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

10 个评论
显示 8更早的评论隐藏 8更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Error using pdepe function in matlab.

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

10 个评论 显示 8更早的评论隐藏 8更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论