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

3 views (last 30 days)
Han F on 20 Nov 2021
Commented: Han F on 23 Nov 2021
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);
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
##### 2 CommentsShowHide 1 older comment
Han F on 22 Nov 2021
Hi Bill,
q and qH are functions of u1 and u2 (q(u1,u2), qH(u1,u2)). There are 2 independent variables (u1, u2).

Bill Greene on 22 Nov 2021
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:
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 on 22 Nov 2021
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

Bill Greene on 23 Nov 2021
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.
##### 2 CommentsShowHide 1 older comment
Han F on 23 Nov 2021
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?

### Categories

Find more on PDE Solvers in Help Center and File Exchange

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!