Boundary using system of pdepes

I have the following pdes systèmes wirh four membres:R_j*dC_j/dt = 1000*(1+0.0025*x)^2 - 100*(1+0.0025*x) - R_j*C_j + R_(j-1)*C_(j-1) with initial and boundary conditions. C_j(x,t=0)=0;. 100*C_j(x=0,t)-1000*C_j(x=0,t)=f_j(t) dC_j(x=L,t)=0 R_j =[10000, 14000, 50000, 500]. I try to solve it with pdepe solver but i obtained results differents to another method. I would like to know if it is because of the definition of boundary conditions in m'y pdepe code. Heure ils the pdepe code

function pdex4 m = 0; x = linspace(0,250,100); t = linspace(0,10000,100);

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3 = sol(:,:,3); u4 = sol(:,:,4);

figure surfc(x,t,u1) title('u1(x,t)') xlabel('Distance x') ylabel('Time t')

figure surfc(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t')

figure surfc(x,t,u3) title('u3(x,t)') xlabel('Distance x') ylabel('Time t')

figure surfc(x,t,u4) title('u4(x,t)') xlabel('Distance x') ylabel('Time t')

figure plot(x,u1(10,:), x,u2(10,:), x,u3(10,:), x,u4(10,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)') % -------------------------------------------------------------- function [c,f,s] = pdex4pde(x,t,u,DuDx) c = [10000; 14000; 50000; 500]; f = [1000*((1+0.0025*x)^2); 1000*((1+0.0025*x)^2); 1000*((1+0.0025*x)^2); 1000*((1+0.0025*x)^2)] .* DuDx ... -[100*(1+0.0025*x); 100*(1+0.0025*x); 100*(1+0.0025*x); 100*(1+0.0025*x)].* u; s = [-0.0079*10000*u(1); 0.0079*10000*u(1)-0.0000028*14000*u(2); 0.0000028*14000*u(2)-0.0000087*50000*u(3); 0.0000087*50000*u(3)-0.00043*500*u(4)]; % -------------------------------------------------------------- function u0 = pdex4ic(x); u0 = [0; 0; 0; 0]; % -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) pl = [1.25*exp(-0.0089*t); 1.25044*(exp(-0.0010028*t)-exp(-0.0089*t)); (0.443684*1e-3*exp(-0.0089*t)+0.593431*exp(-0.0010028*t)-0.593874*exp(-0.0010087*t)); (-0.51674*1e-6*exp(-0.0089*t)+0.120853*1e-1*exp(-0.0010028*t)-0.122637*1e-1*exp(-0.0010087*t)+0.178925*1e-3*exp(-0.00143*t))]; ql = [1; 1; 1; 1]; pr = [100*(1+0.0025*100)*ur(1); 100*(1+0.0025*100)*ur(2); 100*(1+0.0025*100)*ur(3); 100*(1+0.0025*100)*ur(4)]; qr = [1; 1; 1; 1];

[Torsten]

Please format your code - it's unreadable.

And as noted before, "pdepe" is not suited to solve PDEs that don't contain 2nd order derivatives in space.

[Thomas TJOCK-MBAGA]

Sorry i have just forgotten the second order derivative.

Here attached the pdf file of the pdepes code The pde IS R_j*dC_j/dt=d/dx[1000*(1+0.0025)^2*dC_j/dx- 100*(1+0.0025*x)] - R_j*C_j + R{j-1}C{j-1}

[Torsten]

If you include MATLAB code, then it is helpful if you supply the .m file so that participants are able to execute it.

In your case, it will also be necessary to supply the mathematical problem formulation in order to see whether there are errors in the formulation as MATLAB code.

[Thomas TJOCK-MBAGA]

Here attached the problem and the MATLAB code

[Thomas TJOCK-MBAGA]

Here attached the .m file

m = 0;
x = linspace(0,250,100);
t = linspace(0,10000,100);

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);

figure
surfc(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surfc(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surfc(x,t,u3)
title('u3(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surfc(x,t,u4)
title('u4(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
plot(x,u1(10,:), x,u2(10,:), x,u3(10,:), x,u4(10,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')

% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [10000; 14000; 50000; 500]; 
f = [1000*((1+0.0025*x)^2); 1000*((1+0.0025*x)^2); 1000*((1+0.0025*x)^2); 1000*((1+0.0025*x)^2)] .* DuDx ...
    -[100*(1+0.0025*x); 100*(1+0.0025*x); 100*(1+0.0025*x); 100*(1+0.0025*x)].* u; 
s = [-0.0079*10000*u(1); 0.0079*10000*u(1)-0.0000028*14000*u(2); 0.0000028*14000*u(2)-0.0000087*50000*u(3); 0.0000087*50000*u(3)-0.00043*500*u(4)]; 
end
% --------------------------------------------------------------
function u0 = pdex4ic(x);
u0 = [0; 0; 0; 0]; 
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [1.25*exp(-0.0089*t); 1.25044*(exp(-0.0010028*t)-exp(-0.0089*t)); (0.443684*1e-3*exp(-0.0089*t)+0.593431*exp(-0.0010028*t)-0.593874*exp(-0.0010087*t)); (-0.51674*1e-6*exp(-0.0089*t)+0.120853*1e-1*exp(-0.0010028*t)-0.122637*1e-1*exp(-0.0010087*t)+0.178925*1e-3*exp(-0.00143*t))]; 
ql = [1; 1; 1; 1]; 
pr = [100*(1+0.0025*100)*ur(1); 100*(1+0.0025*100)*ur(2); 100*(1+0.0025*100)*ur(3); 100*(1+0.0025*100)*ur(4)]; 
qr = [1; 1; 1; 1]; 
end

[Torsten]

Your PDEs in the mathematical script are incorrect. The convective flux must be u*C, not simply u. But in your code, they seem to be correct.

Which boundary conditions did you set in the "different method" you use ?

[Thomas TJOCK-MBAGA]

I use the same boundary conditions as in the mathematical formulation. I would like to know if the conditions used in the pdepes script deal with the boundary conditions in the mathematical formulation if the problème.

[Torsten]

The boundary conditions in your pdepe code are

at x = 0 (substitute x = 0 in the equations, and the gradients are evaluated at x = 0):

1.25*exp(-0.0089*t) + 1000*((1+0.0025*x)^2)*dC1/dx - 100*(1+0.0025*x)*C1 = 0
1.25044*(exp(-0.0010028*t)-exp(-0.0089*t)) + 1000*((1+0.0025*x)^2)*dC2/dx - 100*(1+0.0025*x)*C2 = 0
(0.443684*1e-3*exp(-0.0089*t)+0.593431*exp(-0.0010028*t)-0.593874*exp(-0.0010087*t)) + 1000*((1+0.0025*x)^2)*dC3/dx -  100*(1+0.0025*x)*C3 = 0
(-0.51674*1e-6*exp(-0.0089*t)+0.120853*1e-1*exp(-0.0010028*t)-0.122637*1e-1*exp(-0.0010087*t)+0.178925*1e-3*exp(-0.00143*t)) + 1000*((1+0.0025*x)^2)*dC4/dx - 100*(1+0.0025*x)*C4 = 0

at x = 250 (substitute x = 250 in the equations, and the gradients are evaluated at x = 250):

100*(1+0.0025*100)*C1 + 1000*((1+0.0025*x)^2)*dC1/dx - 100*(1+0.0025*x)*C1 = 0
100*(1+0.0025*100)*C2 + 1000*((1+0.0025*x)^2)*dC2/dx - 100*(1+0.0025*x)*C2 = 0
100*(1+0.0025*100)*C3 + 1000*((1+0.0025*x)^2)*dC3/dx -  100*(1+0.0025*x)*C3 = 0
100*(1+0.0025*100)*C4 + 1000*((1+0.0025*x)^2)*dC4/dx - 100*(1+0.0025*x)*C4 = 0

[Thomas TJOCK-MBAGA]

Ok, but to obtained these what values for P1, q1, pr and qr am i supposed to used ???

[Torsten]

I don't understand what you mean. Where do P1, q1, pr or qr appear in the 8 equations I wrote ? I just wrote out

pl + ql*f  = 0

and

pr + qr*f = 0

for the values of pl, ql, pr and qr you set in your code.

The general boundary condition for pdepe reads

for the left boundary point:

pl + ql*f = 0

for the right boundary point

pr + qr*f = 0

where f is the flux you defined in "pdex4pde".

You are given xl, xr, ul, ur and t, and you are supposed to return pl, ql, pr, qr to define the boundary conditions of the functions at the left and right boundary points.

If you still have problems, tell me a condition you want to set and I will try to derive p and q from your given f.

[Thomas TJOCK-MBAGA]

Here are the conditions -1000*dC_j(x=0,t)/dt+100*C_j(x=0,t) = f_j(t) dC_j(x=L=250,t)=0 With f_j (t) given in the code

[Torsten]

-1000*dC_j(x=0,t)/dt+100*C_j(x=0,t) = f_j(t)

Are you sure the first derivative is with respect to t ?

If you write your first condition (f from "pdex4pde") as

-f = f_j(t)

or

f + f_j(t) = 0

you can see that you have to set

ql = [1 1 1 1].'
pl = [f_1(t) f_2(t) f_3(t) f_4(t)].'

For the second it's more difficult.

pr + qr*(1000*((1+0.0025*xr)^2) * dC/dx - 100*(1+0.0025*xr) * C) 

should become

dC/dx.

So choose

qr = 1/(1000*((1+0.0025*xr)^2))*[1 1 1 1].';
pr = 100*(1+0.0025*xr) * ur / (1000*((1+0.0025*xr)^2))

[Thomas TJOCK-MBAGA]

Thanks you

[Torsten]

I just noticed that the normalization by 1000*((1+0.0025*xr)^2) was not necessary because

1000*((1+0.0025*xr)^2) * dC/dx = 0

implies

dC/dx = 0

Thus your boundary conditions in the code are in correspondence with your mathematical equations if you change the 100 by xr ( = 250) in

pr = [100*(1+0.0025*100)*ur(1); 100*(1+0.0025*100)*ur(2); 100*(1+0.0025*100)*ur(3); 100*(1+0.0025*100)*ur(4)];