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

1 view (last 30 days)
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
  2 Comments
Bill Greene
Bill Greene on 21 Nov 2021
What are q and qH? It looks like you have four dependent variables and only two equations.
Han F
Han F on 22 Nov 2021
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).

Sign in to comment.

Accepted Answer

Bill Greene
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:
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
  1 Comment
Han F
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

Sign in to comment.

More Answers (1)

Bill Greene
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 Comments
Han F
Han F on 23 Nov 2021
Thank you, Bill. I really appreciate the idea you provided. Especially using inverse to rewrite equation coefficients, genius.
Han F
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?

Sign in to comment.

Tags

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!