Error using Time dependant Boundary Conditions with PDEPE

I have defined a system of six pde equations to solve a diffussion model. As boundary condition on the right side I need to introduce experimental values that are differerent for each time and cannot be adjusted to a curve. Thus, they have to be introduce as a vector of numbers.... Is this possible??

My code for the Boundary Conditions is the following:

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions

global data Nt

pl = [0; 0; 0; 0; 0; 0];

ql = [1; 1; 1; 1; 1; 1];

pr = [data(1,t)-ur(1); data(2,t)-ur(2); data(3,t)-ur(3); data(4,t)-ur(4); (data(2,t)+data(3,t)+2*data(4,t))-ur(5); (data(1,1)-data(1,t)-data(2,t)-data(3,t)-data(4,t))-ur(6)];

qr = [0; 0; 0; 0; 0; 0];

end

In the previous code, the matrix "data" contains the experimental data of 4 of the dependant variables. The BC for the last two dependant variables in the term "pr" can be calculated from the other.

Any help would be great. Thanks!

Answer 1

[Torsten]

If "time" is the array corresponding to the data matrix, use

data_actual = interp1(time,data,t) 

and set

pr = [data_actual(1)-ur(1); ...
      data_actual(2)-ur(2); ...
      data_actual(3)-ur(3); ...
      data_actual(4)-ur(4); ...
      (data_actual(2)+data_actual(3)+2*data_actual(4))-ur(5); ...
      (data(1,1)-data_actual(1)-data_actual(2)-data_actual(3)-data_actual(4))-ur(6)];

[Pablo Gutiérrez Sánchez]

Hi Torsten,

First of all, thank you for your time and your quick reply!

Regarding your suggestion, I think I am doing something wrong, because it is not working. The error message is the following:

Error using interp1>reshapeAndSortXandV

LENGTH(X) and SIZE(V,1) must be the same.

Error in interp1 (line 128)

[X,V,orig_size_v] = reshapeAndSortXandV(X,V);

Error in soloecuacion>pdebc (line 288)

data_actual = interp1(timedata,data,t)

Error in pdepe (line 250)

[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});

Error in soloecuacion (line 59)

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

I think it is easier if I share all the code:

close all
clear all
global params0 timedata data h2o_exp2 hums_exp2

L = 1;     


timedata = [0 15 30 45 60 75 105 135 195 255 315 375];
t = [0 15 30 45 60 75 105 135 195 255 315 375];
Nt=length(t);
Nx=length(t);
x = linspace(0, L, Nx);

% Experimental data

stol_exp2=[0.013500576, 0.009736718, 0.007708470, 0.006096486, 0.004548828, 0.003667008, 0.002276970, 0.001318896, 0.000387234, 0.000129504, 0.000017466, 0.000002769];

stanG_exp2=[0.000000000, 0.000219786, 0.000885996, 0.001949614, 0.002742847, 0.003480738, 0.004303486, 0.005244297, 0.005745536, 0.006246775, 0.006380650, 0.006254681]; %viavle

stanB_exp2=[0.000000000, 0.000000000, 0.000046382, 0.000158647, 0.000257208, 0.000338376, 0.000461182, 0.000618774, 0.000785469, 0.000834647, 0.000857612, 0.000861424]; %no viable

ibda_exp2=[0.000000000, 0.000256459, 0.000549803, 0.000671379, 0.000838870, 0.001061570, 0.001234599, 0.001528664, 0.002311477, 0.002547755, 0.003142632, 0.003204195];



hums_exp2=0;
h2o_exp2=0.416666667;


data = [stol_exp2; stanG_exp2; stanB_exp2; ibda_exp2];


k10=134;
k20=216;
k30=148;
k40=280;
Kstol0=200; 
KstanV0=263; 
KstanNV0=208; 
Kiso0=21; 
Kh2o0=250; 

params0=[k10 k20 k30 k40 Kstol0 KstanV0 KstanNV0 Kiso0 Kh2o0]; 



options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7) 

options = struct with fields:
              AbsTol: 1.0000e-04
                 BDF: []
              Events: []
         InitialStep: 1.0000e-07
            Jacobian: []
           JConstant: []
            JPattern: []
                Mass: []
        MassSingular: []
            MaxOrder: []
             MaxStep: []
         NonNegative: []
         NormControl: 'off'
           OutputFcn: []
           OutputSel: []
              Refine: []
              RelTol: 1.0000e-04
               Stats: []
          Vectorized: []
    MStateDependence: []
           MvPattern: []
        InitialSlope: []

m=2

m = 2

sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t,options);
u1 = sol(:,:,1); %stol
u2 = sol(:,:,2); 
u3 = sol(:,:,3); 
u4 = sol(:,:,4); 
u5 = sol(:,:,5); 
u6 = sol(:,:,6); 
x;
t;

%compuestos clave

u1xinicial=[t;u1(:,1).'];
u1xfin=[t;u1(:,Nx).'];

u1tinicial=[x;u1(1,:)];
u1tfinal=[x;u1(Nt,:)];

u2xinicial=[t;u2(:,1).'];
u2xfin=[t;u2(:,Nx).'];

u2tinicial=[x;u2(1,:)];
u2tfinal=[x;u2(Nt,:)];

u3xinicial=[t;u3(:,1).'];
u3xfin=[t;u3(:,Nx).'];

u3tinicial=[x;u3(1,:)];
u3tfinal=[x;u3(Nt,:)];

u4xinicial=[t;u4(:,1).'];
u4xfin=[t;u4(:,Nx).'];

u4tinicial=[x;u4(1,:)];
u4tfinal=[x;u4(Nt,:)];




% graphs

figure (1);
subplot(1,6,1);
surf(x, t, u1);
xlabel('x');
ylabel('t');
zlabel('u1');
title('Variable u1');
subplot(1,6,2);
surf(x, t, u2);
xlabel('x');
ylabel('t');
zlabel('u2');
title('Variable u2');
subplot(1,6,3);
surf(x, t, u3);
xlabel('x');
ylabel('t');
zlabel('u3');
title('Variable u3');
subplot(1,6,4);
surf(x, t, u4);
xlabel('x');
ylabel('t');
zlabel('u4');
title('Variable u4');
subplot(1,6,5);
surf(x, t, u5);
xlabel('x');
ylabel('t');
zlabel('u3');
title('Variable u5');
subplot(1,6,6);
surf(x, t, u6);
xlabel('x');
ylabel('t');
zlabel('u6');
title('Variable u6');

figure (2);
subplot(1,6,1);
scatter(t, data (1,:));
hold on
plot(t, u1xfin(2,:))
xlabel('t');
ylabel('u1');
title('Datos u1 en X fin');
subplot(1,6,2);
scatter(t, data (2,:));
hold on
plot(t, u2xfin(2,:))
xlabel('t');
ylabel('u2');
title('Datos u2 en X fin');
subplot(1,6,3);
scatter(t, data (3,:));
hold on
plot(t, u3xfin(2,:))
xlabel('t');
ylabel('u3');
title('Datos u3 en X fin');
subplot(1,6,4);
scatter(t, data (4,:));
hold on
plot(t, u4xfin(2,:))
xlabel('t');
ylabel('u4');
title('Datos u4 en X fin');
subplot(1,6,5);
plot(t, u5xfin(2,:))
xlabel('t');
ylabel('u5');
title('Datos u5 en X fin');
subplot(1,6,6);
plot(t, u6xfin(2,:))
xlabel('t');
ylabel('u6');
title('Datos u6 en X fin');

figure (3);
subplot(1,6,1);
plot(x, u1tfinal(2,:))
xlabel('x');
ylabel('u1');
title('Datos u1 en t fin');
subplot(1,6,2);
plot(x, u2tfinal(2,:))
xlabel('x');
ylabel('u2');
title('Datos u2 en t fin');
subplot(1,6,3);
plot(x, u3tfinal(2,:))
xlabel('x');
ylabel('u3');
title('Datos u3 en t fin');
subplot(1,6,4);
plot(x, u4tfinal(2,:))
xlabel('x');
ylabel('u4');
title('Datos u4 en t fin');
subplot(1,6,5);
plot(x, u5tfinal(2,:))
xlabel('x');
ylabel('u5');
title('Datos u5 en t fin');
subplot(1,6,6);
plot(x, u6tfinal(2,:))
xlabel('x');
ylabel('u6');
title('Datos u6 en t fin');

function [c,f,s] = pdefun(x,t,u,dudx,params0) 
global params0

DstanV=(4.244*10^-5)*60;
DstanNV=(4.244*10^-5)*60;
Rp=0.05; 
cat=450; 

%Parameters

k1=params0(1);
k2=params0(2);
k3=params0(3);
k4=params0(4);
Kstol=params0(5); 
KstanV=params0(6); 
KstanNV=params0(7); 
Kiso=params0(8); 
Kh2o=params0(9); 
 
% kinetics

r1 = (k1*Kstol*u(1))/((1+KstanNV*u(3)+Kiso*u(4))^2); 
r2 = (k2*Kstol*u(1))/((1+Kstol*u(1)+KstanNV*u(3)+Kiso*u(4))^2); 
r3 = (k3*KstanV*u(2))/((1+KstanNV*u(3)+Kiso*u(4))^2); 
r4 = (k4*Kstol*u(1))/((1+Kstol*u(1)+Kiso*u(4))^2); 



Rstol=(r2-r4)*cat;
RstanV=(+r3)*cat;
RstanNV=r2*cat;
Ribda=r3*cat;
Ragua=(r4+r3)*cat;
Rhums=r4*cat;

%pdepe

c = [1; 1; 1; 1; 1; 1];
f = [0; DstanV/Rp^2; DstanNV/Rp^2; 0; 0; 0] .* dudx;
s = [Rstol; RstanV; RstanNV; Ribda; Ragua; Rhums];

end


% ---------------------------------------------
function u0 = pdeic(x) % Initial Conditions
global data h2o_exp2 hums_exp2

    u0 = [data(1,1); data(2,1); data(3,1); data(4,1); h2o_exp2; hums_exp2];
end

% ---------------------------------------------
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
global timedata data h2o_exp2 hums_exp2

pl = [0; 0; 0; 0; 0; 0];
ql = [1; 1; 1; 1; 1; 1];
data_actual = interp1(timedata.',data.',t);
pr = [data_actual(1)-ur(1); ...
      data_actual(2)-ur(2); ...
      data_actual(3)-ur(3); ...
      data_actual(4)-ur(4); ...
      (data_actual(2)+data_actual(3)+2*data_actual(4))-ur(5); ...
      (data(1,1)-data_actual(1)-data_actual(2)-data_actual(3)-data_actual(4))-ur(6)];
%pr = [data(1,t)-ur(1); data(2,t)-ur(2); data(3,t)-ur(3); data(4,t)-ur(4); data(2,t)+data(3,t)+2*data(4,t)+h2o_exp2-ur(5); data(1,1)-data(1,t)-data(2,t)-data(3,t)-data(4,t)+humins_exp2-ur(6)];
qr = [0; 0; 0; 0; 0; 0];
end

As you can see, the boundary conditions for u1 are time dependant values and are included in data(1,:); for u2 are included in data(2,:)... An the time vector is specified as "timedata" or "t" (both are the same and defined at the beginning of the code).

Now you have the whole code, do you have any clue of the error?

As said before, thank you!!

[Torsten]

Why didn't you include my suggested changes in your code ? It works now (at least technically) (see above).

[Pablo Gutiérrez Sánchez]

Thank your Torsten. But I still have some doubts... As you can see in the code, the term "f" in @pdefun specifies that u1, u4, u5 and u6 do not change with the independat variable "x" (f=0 in all of them), but Figure 3 seems not to follow this behaviour. That is, Figure 3 is correct for the two dependant variables u2 and u3, but the rest shoud show a constant value along x axis (and this value should be equal to the boundary condition at the right side). At this point, you will have realized that I am trying to solve a system of two PDEs (u2 and u3) and four ODEs (u1, u4,u5,u6). The thing is that for the last group I only need to specify the initial condition, but I think the code is using all the experimental points. How could I fix this issue?

[Torsten]

As you can see in the code, the term "f" in @pdefun specifies that u1, u4, u5 and u6 do not change with the independat variable "x" (f=0 in all of them), but Figure 3 seems not to follow this behaviour.

The array "s" for u1, u4, u5 and u6 varies with x because u2 and u3 do. So I don't understand why you expect a constant profile for u1, u4, u5 and u6.

And strictly speaking, using "pdepe" for f = 0 is not allowed. "pdepe" solves parabolic-elliptic PDEs, not ODEs. You should at least specify f (i) = (small_value)*dudx(i) for i = 1, 4, 5 and 6.

[Pablo Gutiérrez Sánchez]

Many thanks for your help Torsten!!