How to deal these type of errors from pdepe? "Error in pdepe/pdeodes (line 359)"

Hello,

I am trying to solve a system of pdes with the following form:

I wonder if it is possible to solve this system with pdepe function?

Currently, I got errors like:

Unable to perform assignment because the size of the left side is 26-by-1 and the size of the right side is 26-by-26.

Error in pdepe/pdeodes (line 359)

up(:,ii) = ((xim(ii) * fR - xim(ii-1) * fL) + ...

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 152)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in pdepe (line 289)

[t,y] = ode15s(@pdeodes,t,y0(:),opts);

Error in PDEPE_solver (line 8)

sol = pdepe(m,@pdefun,@pdeic,@pdebc,z,t);

Many thanks in advance!

Han

clear all
global zz
z = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
zz=z;
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,z,t);

function [c,f,s] = pdefun(z,t,u,dudz) %(S is u(1), T is u(2))
global zz
[~,n]=size(zz);



[M, R, rho, lamta, he, g, yeeta, ksat, thetasat, thetares, thetan,thetao, rhob, thetac,Cv,CW]=constant(n);

theta=(thetasat-thetares).*u(1)+thetares;

thao=((thetasat-theta).^(7/3))./((thetasat).^2);

Dv=2.12e-5.*10.^5./101325.*((u(2)+273.15)./273.15).^1.88.*thao.*(thetasat-theta);

kH=(0.65-0.78.*rhob./1000+0.6.*(rhob./1000).^2)+(1.06.*(rhob./1000)).*theta-...
        ((0.65-0.78.*rhob./1000+0.6.*(rhob./1000).^2)-(0.03+0.1*(rhob./1000)^2)).*exp(-((1+2.6*(thetac)^(-0.5)).*theta).^(4));

Csoil=(1.92*thetan+2.51*thetao+4.18*theta)*10^6;

lamtaE=(2.501-2.361*10^(-3)*15)*10^6;

phaie=ksat.*he./(1-lamta.*yeeta);

cvsat=610.78.*exp(17.27.*u(2)./(u(2)+237.3)).*M./R./rho./(u(2)+273.15);

h=(u(1)).^(-1./lamta).*he;

hr=exp(h.*g.*M./R./(u(2)+273.15));

dhrdS=-((exp((g.*he.*M.*u(1).^(-1./lamta))./(R.*(273.15 + u(2)))).*g.*he.*M.*u(1).^(-1 - 1./lamta))./(lamta.*R.*(273.15+u(2))));

dhrdT=-((exp((g.*he.*M.* u(1).^(-1./lamta))./(R.*(273.15 + u(2)))).*g.*he.*M.*u(1).^(-1./lamta))./(R.* (273.15 + u(2)).^2));

dcvsatdT=-(610.78.*exp((17.27.*u(2))./(237.3+u(2))).*M.*(-1.0631.*10^6-3623.57.*u(2)+u(2).^2))./(R.*rho.*(237.3+u(2)).^2.*(273.15+u(2)).^2);

dphaidS=(phaie.*(lamta.*yeeta - 1))./(u(1).^(1./lamta + 1).*lamta.*(1./u(1).^(1./lamta)).^(lamta.*yeeta));

k=ksat.*(u(1)).^yeeta;

dkdS=u(1).^(yeeta - 1).*ksat.*yeeta;

%-----------

M1=(thetasat-thetares).*(1-cvsat.*hr)+(thetasat-thetares).*(1-u(1)).*cvsat.*dhrdS;

M2=(thetasat-thetares).*(1-u(1)).*dcvsatdT.*hr+(thetasat-thetares).*(1-u(1)).*cvsat.*dhrdT;

M3=(thetasat-thetares).*rho.*lamtaE.*(1-u(1)).*cvsat.*dhrdS-(thetasat-thetares).*rho.*lamtaE.*cvsat.*dhrdT;

M4=Csoil+(thetasat-thetares).*rho.*lamtaE.*(1-u(1)).*dcvsatdT.*hr+(thetasat-thetares).*rho.*lamtaE.*(1-u(1)).*cvsat.*dhrdT;

%--------

M11=-(-dphaidS-Dv.*cvsat.*dhrdS);

M22=Dv.*hr.*dcvsatdT+Dv.*cvsat.*dhrdT;

M33=-(-dphaidS.*CW.*u(2)-(Cv.*u(2)+rho.*lamtaE).*Dv.*cvsat.*dhrdS);

M44=-(-kH-(Cv.*u(2)+rho.*lamtaE).*Dv.*hr.*dcvsatdT-(Cv.*u(2)+rho.*lamtaE).*Dv.*cvsat.*dhrdT);

%-----------

M111=-dkdS;

M222(1:n,1)=0;

M333=-CW.*u(2).*dkdS;

M444=-CW.*k;

%-------------
c_index_row(1:4:4*n-3,1)=1:2:2*n-1;
c_index_row(2:4:4*n-2,1)=2:2:2*n;
c_index_row(3:4:4*n-1,1)=1:2:2*n-1;
c_index_row(4:4:4*n,1)=2:2:2*n;

c_index_column(1:2:4*n-1,1)=1:2*n;
c_index_column(2:2:4*n,1)=1:2*n;

c_index_M=zeros(4*n,1);
c_index_M(1:4:4*n-3)=M1;
c_index_M(3:4:4*n-1)=M2;
c_index_M(2:4:4*n-2)=M3;
c_index_M(4:4:4*n)=M4;

c=sparse(c_index_row,c_index_column,c_index_M);


%-------------------f matrix


f_index_M=zeros(4*n,1);
f_index_M(1:4:4*n-3)=M11;
f_index_M(3:4:4*n-1)=M22;
f_index_M(2:4:4*n-2)=M33;
f_index_M(4:4:4*n)=M44;

f = sparse(c_index_row,c_index_column,f_index_M)* dudz;
%-------------------s matrix


s_index_M=zeros(4*n,1);
s_index_M(1:4:4*n-3)=M111;
s_index_M(3:4:4*n-1)=M222;
s_index_M(2:4:4*n-2)=M333;
s_index_M(4:4:4*n)=M444;

s = sparse(c_index_row,c_index_column,s_index_M);
end

function u0 = pdeic(z)
global zz
[~,n]=size(zz);
S_initial(1:n,1)=0.2;
T_initial(1:n,1)=10;

u0(1:2:2*n-1,1) = S_initial;%[S_initial; T_initial];
u0(2:2:2*n,1) = T_initial;
end

function [pl,ql,pr,qr] = pdebc(zl,ul,zr,ur,t)
global zz

[~,n]=size(zz);

S0=0.2; T0=10; 
Sn=0.1; Tn=10;
upper_boundary(1:2:2*n-1,1)=S0;
upper_boundary(2:2:2*n,1)=T0;

lower_boundary(1:2:2*n-1,1)=Sn;
lower_boundary(2:2:2*n,1)=Tn;


pl = [ul-upper_boundary];
ql = zeros(2*n,1);%[0; 0];
pr = [ur-lower_boundary];
qr = zeros(2*n,1);%[0; 0];
end

function [M, R, rho, lamta, he, g, yeeta, ksat, thetasat, thetares, thetan,thetao, rhob, thetac,Cv,CW]=constant(n)
M=0.018;
R=8.316;
rho=1000;
lamta(1:n,1)=0.15;
he(1:n,1)=-1/3.07;
g=9.8;
yeeta(1:n,1)=2/lamta+3;
ksat(1:n,1)=1.16e-5;
thetasat(1:n,1)=0.48;
thetares(1:n,1)=0.04;

rhob=1300;

Cv=1.8*10^6;                                                               % volumetric heat capacity of vapor, J/m3/K
CW=4.18*10^6;

thetac=0.2;
thetam=0.393;
thetaq=0.243;
thetan=0.529;
thetao=0;

end

[Bill Greene]

What are q and qH? It looks like you have four dependent variables and only two equations.

[Han F]

Hi Bill,

Thank your for your reply.

q and qH are functions of u1 and u2 (q(u1,u2), qH(u1,u2)). There are 2 independent variables (u1, u2).

Answer 1

[Bill Greene]

pdepe does not accept a non-diagonal mass matrix. But often you can deal with this by calculating the inverse of the mass matrix and using that to transform the equations. Here is a simple two-equation example:

function matlabAnswers_11_22_2021

L=1;

n=30;

n2 = ceil(n/2);

x = linspace(0,L,n);

tf=.1;

t = linspace(0,tf,30);

M=[11 3; 3 2];

sola=diffusionEquationCoupledMassAnalSoln(t, x, M, @icFunc);

pdef = @(x,t,u,DuDx) pdeFunc(x,t,u,DuDx,M);

icf = @(x) icFunc(x, L);

bcFunc = @(xl,ul,xr,ur,t) bcDirichlet(xl,ul,xr,ur,t);

m=0;

opts=struct;

opts.RelTol=1e-5;

opts.AbsTol=1e-7;

sol = pdepe(m, pdef,icf,bcFunc,x,t,opts);

u1=sol(:,:,1);

u2=sol(:,:,2);

u1a=sola(:,:,1);

u2a=sola(:,:,2);

figure; plot(t, u1(:, n2), t, u2(:, n2), ...

t, u1a(:,n2), '+', t, u2a(:,n2), 'o'); grid on;

xlabel 'time'

title 'solution at center';

legend('u1', 'u2', 'u1,anal', 'u2,anal');

figure; plot(x, u1(end, :), x, u2(end, :), ...

x, u1a(end, :), '+', x, u2a(end, :), 'o'); grid on;

xlabel 'x'

title('solution at final time');

legend('u1', 'u2', 'u1,anal', 'u2,anal');

err=max(abs(sol(:)-sola(:)));

fprintf('Solution: Time=%g, u_center=%g, v_center=%g, err=%10.2e\n', ...

t(end), u1(end,n2), u2(end,n2), err);

end

function [c,f,s] = pdeFunc(x,t,u,DuDx,M)

nx=length(x);

c = ones(2,nx);

f = inv(M)*DuDx;

s=zeros(2,nx);

end

function u0 = icFunc(x,L)

u0 = [sin(pi*x/L); sin(pi*x/L)];

end

% --------------------------------------------------------------

function [pl,ql,pr,qr] = bcDirichlet(xl,ul,xr,ur,t)

pl = ul;

ql = [0 0]';

pr = ur;

qr = [0 0]';

end

function u=diffusionEquationCoupledMassAnalSoln(t, x, M, icFunc)

nt=length(t);

nx=length(x);

nm=size(M,1);

iM=inv(M);

[vec,val]=eig(iM);

%iM-vec*val*vec'

L=x(end);

nInt=1000;

xInt = linspace(0,L,nInt);

u0=icFunc(xInt,L);

u0=vec'*u0; % transform IC to modal space

D1=2/L*trapz(xInt/L, (u0.*sin(pi*xInt/L))');

u=zeros(nt,nx,nm);

w=zeros(nt,nx,nm);

for i=1:nm % solve the uncopupled equations

et=exp(-pi^2*val(i,i)*t(:)/L^2);

v=(D1(i)*sin(pi*x/L));

w(:,:,i) = et*v;

end

for i=1:nt

u(i,:,:) = (vec*squeeze(w(i,:,:))')';

end

end

[Han F]

Thank you Bill, I can get solutions now. However, the M matrix in your example is constant for whole x, I wonder how to deal with non-constant M matrix.

I mean, if I divide whole profile to 13 layers, coefficients for each layer are different. Finally, I want get the size of sol is length(t)*length(x)*2 (u1 and u2 profile at each time step). Currently, sol is length(t)*length(x)*26, it seems each equations is solved at whole x.

I attach my code here, I appreciate your help very much.

clear all
global zz
z = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
zz=z;
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,z,t);

function [c,f,s] = pdefun(z,t,u,dudz) %(S is u(1), T is u(2))
global zz
[~,n]=size(zz);



[M, R, rho, lamta, he, g, yeeta, ksat, thetasat, thetares, thetan,thetao, rhob, thetac,Cv,CW]=constant(n);

theta=(thetasat-thetares).*u(1)+thetares;

thao=((thetasat-theta).^(7/3))./((thetasat).^2);

Dv=2.12e-5.*10.^5./101325.*((u(2)+273.15)./273.15).^1.88.*thao.*(thetasat-theta);

kH=(0.65-0.78.*rhob./1000+0.6.*(rhob./1000).^2)+(1.06.*(rhob./1000)).*theta-...
        ((0.65-0.78.*rhob./1000+0.6.*(rhob./1000).^2)-(0.03+0.1*(rhob./1000)^2)).*exp(-((1+2.6*(thetac)^(-0.5)).*theta).^(4));

Csoil=(1.92*thetan+2.51*thetao+4.18*theta)*10^6;

lamtaE=(2.501-2.361*10^(-3)*15)*10^6;

phaie=ksat.*he./(1-lamta.*yeeta);

cvsat=610.78.*exp(17.27.*u(2)./(u(2)+237.3)).*M./R./rho./(u(2)+273.15);

h=(u(1)).^(-1./lamta).*he;

hr=exp(h.*g.*M./R./(u(2)+273.15));

dhrdS=-((exp((g.*he.*M.*u(1).^(-1./lamta))./(R.*(273.15 + u(2)))).*g.*he.*M.*u(1).^(-1 - 1./lamta))./(lamta.*R.*(273.15+u(2))));

dhrdT=-((exp((g.*he.*M.* u(1).^(-1./lamta))./(R.*(273.15 + u(2)))).*g.*he.*M.*u(1).^(-1./lamta))./(R.* (273.15 + u(2)).^2));

dcvsatdT=-(610.78.*exp((17.27.*u(2))./(237.3+u(2))).*M.*(-1.0631.*10^6-3623.57.*u(2)+u(2).^2))./(R.*rho.*(237.3+u(2)).^2.*(273.15+u(2)).^2);

dphaidS=(phaie.*(lamta.*yeeta - 1))./(u(1).^(1./lamta + 1).*lamta.*(1./u(1).^(1./lamta)).^(lamta.*yeeta));

k=ksat.*(u(1)).^yeeta;

dkdS=u(1).^(yeeta - 1).*ksat.*yeeta;

%-----------

M1=(thetasat-thetares).*(1-cvsat.*hr)+(thetasat-thetares).*(1-u(1)).*cvsat.*dhrdS;

M2=(thetasat-thetares).*(1-u(1)).*dcvsatdT.*hr+(thetasat-thetares).*(1-u(1)).*cvsat.*dhrdT;

M3=(thetasat-thetares).*rho.*lamtaE.*(1-u(1)).*cvsat.*dhrdS-(thetasat-thetares).*rho.*lamtaE.*cvsat.*dhrdT;

M4=Csoil+(thetasat-thetares).*rho.*lamtaE.*(1-u(1)).*dcvsatdT.*hr+(thetasat-thetares).*rho.*lamtaE.*(1-u(1)).*cvsat.*dhrdT;

%--------

M11=-(-dphaidS-Dv.*cvsat.*dhrdS);

M22=Dv.*hr.*dcvsatdT+Dv.*cvsat.*dhrdT;

M33=-(-dphaidS.*CW.*u(2)-(Cv.*u(2)+rho.*lamtaE).*Dv.*cvsat.*dhrdS);

M44=-(-kH-(Cv.*u(2)+rho.*lamtaE).*Dv.*hr.*dcvsatdT-(Cv.*u(2)+rho.*lamtaE).*Dv.*cvsat.*dhrdT);

%-----------

M111=-dkdS;

M222(1:n,1)=0;

M333=-CW.*u(2).*dkdS;

M444=-CW.*k;

%-------------
c_index_row(1:4:4*n-3,1)=1:2:2*n-1;
c_index_row(2:4:4*n-2,1)=2:2:2*n;
c_index_row(3:4:4*n-1,1)=1:2:2*n-1;
c_index_row(4:4:4*n,1)=2:2:2*n;

c_index_column(1:2:4*n-1,1)=1:2*n;
c_index_column(2:2:4*n,1)=1:2*n;

c_index_M=zeros(4*n,1);
c_index_M(1:4:4*n-3)=M1;
c_index_M(3:4:4*n-1)=M2;
c_index_M(2:4:4*n-2)=M3;
c_index_M(4:4:4*n)=M4;

c_sparse=sparse(c_index_row,c_index_column,c_index_M);
c=diag(ones(2*n,1));
inv_c=c_sparse\c;
c=ones(2*n,1);

%c=c_index_M;

%-------------------f matrix


f_index_M=zeros(4*n,1);
f_index_M(1:4:4*n-3)=M11;
f_index_M(3:4:4*n-1)=M22;
f_index_M(2:4:4*n-2)=M33;
f_index_M(4:4:4*n)=M44;

dudz_M=zeros(4*n,1);
dudz_M(1:4:4*n-3,1)=dudz(1);
dudz_M(2:4:4*n-2,1)=dudz(2);
dudz_M(3:4:4*n-1,1)=dudz(1);
dudz_M(4:4:4*n,1)=dudz(2);

f = sparse(c_index_row,c_index_column,f_index_M);%* dudz_M;
f=inv_c*f*dudz;
%f=f_index_M.*dudz_M;
%-------------------s matrix


s_index_M=zeros(4*n,1);
s_index_M(1:4:4*n-3)=M111;
s_index_M(3:4:4*n-1)=M222;
s_index_M(2:4:4*n-2)=M333;
s_index_M(4:4:4*n)=M444;

u_s=zeros(2*n,1);
u_s(1:2:2*n-1,1)=u(1);
u_s(2:2:n,1)=u(2);
s = sparse(c_index_row,c_index_column,s_index_M);
s =inv_c*s*u_s;
%s= s_index_M;
end

function u0 = pdeic(z)
global zz
[~,n]=size(zz);
S_initial(1:n,1)=0.2;
T_initial(1:n,1)=10;

 u0(1:2:2*n-1,1) = S_initial;%[S_initial; T_initial];
 u0(2:2:2*n,1) = T_initial;

end

function [pl,ql,pr,qr] = pdebc(zl,ul,zr,ur,t)
global zz

[~,n]=size(zz);

S0=0.2; T0=10; 
Sn=0.1; Tn=10;
upper_boundary(1:2:2*n-1,1)=S0;
upper_boundary(2:2:2*n,1)=T0;

lower_boundary(1:2:2*n-1,1)=Sn;
lower_boundary(2:2:2*n,1)=Tn;


pl = ul-upper_boundary;
ql = zeros(2*n,1);%[0; 0];
pr = ur-lower_boundary;
qr = zeros(2*n,1);%[0; 0];
end

function [M, R, rho, lamta, he, g, yeeta, ksat, thetasat, thetares, thetan,thetao, rhob, thetac,Cv,CW]=constant(n)
M=0.018;
R=8.316;
rho=1000;
lamta(1:n,1)=0.15;
he(1:n,1)=-1/3.07;
g=9.8;
yeeta(1:n,1)=2/lamta+3;
ksat(1:n,1)=1.16e-5;
thetasat(1:n,1)=0.48;
thetares(1:n,1)=0.04;

rhob=1300;

Cv=1.8*10^6;                                                               % volumetric heat capacity of vapor, J/m3/K
CW=4.18*10^6;

thetac=0.2;%0.15                                                           % content of clay
thetam=0.393;%0.25; %0.05                                                  % content of other material
thetaq=0.243;%0.25; %0.05                                                  % content of quatze
thetan=0.529;%1-thetasat;                                                  % content of soild
thetao=0;%0.25; %0.05                                                      % content of organic matter% layered_he=                [-1/3.07;-1/2;-1/2.8]; % 3.07,2.1,                      %[-1/3.07;-1/2.07]   [-1/3.07;-1/2.83;-1/3.07];
end

Answer 2

[Bill Greene]

Yes, in my example, the M-matrix was constant but the same idea applies if M is a function of x or the dependent variables. The only problem that would arise is if at some value of x or some values of the dependent variables, the M-matrix becomes singular.

[Han F]

Thank you, Bill. I really appreciate the idea you provided. Especially using inverse to rewrite equation coefficients, genius.

[Han F]

Hi Bill,

A following question is, why sometimes pdepe only return results of first time step? For example, when I set t is linspace(0,1,10), space z= linspace(0,1,100), sometimes size(sol)=1*100*2, while sometimes size(sol)=10*100*2. This relate to initial and boundary conditions?