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

2 次查看（过去 30 天）
Han F
Han F 2021-11-20
评论： Han F 2021-11-23
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 个评论
Bill Greene
Bill Greene 2021-11-21
What are q and qH? It looks like you have four dependent variables and only two equations.
Han F
Han F 2021-11-22
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).

请先登录，再进行评论。

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

采纳的回答

Bill Greene
Bill Greene 2021-11-22
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 个评论
Han F
Han F 2021-11-22
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

请先登录，再进行评论。

更多回答（1 个）

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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by