I am trying to solve the system of coupled partial differential equations described in the attachment using the function pdepe. My code runs into this error:

function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve
c = [1; 1];
f = [-0.6; 1.32].*u;
s = [0; 0];
end
function u0 = pdeic(x) % Initial Conditions
u0 = [sin(pi*x); sin(pi*x)];
end
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
pl = [2.2*ul(2)+ ul(1); 0];
ql = [0; 0];
pr = [0; ur(2)];
qr = [0; 0];
end
x = linspace(0,1,50);
t = linspace(0,2,50);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
surf(x,t,u1)
title('u_(x,t)')
xlabel('Distance x')
ylabel('Time t')

Answer 1

[Torsten]

"pdepe" is a solver for parabolic-elliptic partial differential equations (with a second-order spatial derivative term modelling diffusion in each of the equations).

Your equations are pure hyperbolic in nature (only convection without diffusion) and thus cannot be solved with "pdepe".

The "technical" problems with the solver are that you don't prescribe boundary conditions for the first equation in x = 1 (pr(1) = qr(1) = 0) and for the second equation in x = 0 (pl(2) = ql(2) = 0). That's mathematically correct, but - as said - the solver is designed for second-order PDEs which need boundary conditions at both ends of the integration interval.

Further, f in "pdefun" should be f = [0;0] and s should be s = [-0.6; 1.32].*dudx .

[Yidan]

Thanks for answer!!! I have learnt a lot from your answer.

Sorry, I still have some questions, hope you can help.

1. How to distinguish the convection and diffusion?

2. There is an example in matlab helper that does not have the econd-order spatial derivative term, why the example in the picture works smoothly?

[Torsten]

1. First-order derivative in space is associated with convection, second-order derivative in space is associated with diffusion.
2. There is no MATLAB tool that automatically solves first-order hyperbolic PDEs like "pdepe" for second-order PDEs. You will have to discretize your spatial derivatives and solve the resulting system of ordinary differential equations using ODE15S. Look up "method-of-lines" for more details. Or you can try the "method-of characteristics".

I don't understand the example in the picture. What does "in the" and "body(P t)" mean ?

[Yidan]

Sorry, there are some mistakes in the picture. The right functions are as following:

The link is:求解一维抛物型和椭圆型 PDE - MATLAB pdepe - MathWorks 中国

[Torsten]

The boundary conditions set for this PDE are mathematically incorrect. From the mathematical standpoint, only a condition at x = 1 can be set. The condition at x = 0 is artificial to make "pdepe" work for this problem. I didn't check whether it influences the correct result.

[Torsten]

Here is the code for your problem with the method-of-lines approach:

x = linspace(0,1,500);
t = linspace(0,1,30);
n = numel(x);
M = eye(2*n);
M(1,1) = 0;
M(2*n,2*n) = 0;
y0 = [sin(pi*x),sin(pi*x)];
options = odeset('Mass',M);
[T,Y] = ode15s(@(t,y)fun(t,y,x.'),t,y0,options);
alpha = Y(:,1:n);
beta = Y(:,n+1:2*n);
figure(1)
hold on
plot(x,alpha(end,:),'r')
plot(x,beta(end,:),'b')
hold off
title 'solution at final time'; legend('alpha', 'beta');
grid on

figure(2)
hold on
plot(t,alpha(:,1),'r')
plot(t,beta(:,1),'b')
hold off
title 'solution at left end as a function of time'; legend('alpha', 'beta');
grid on

figure(3)
hold on
plot(t,alpha(:,end),'r')
plot(t,beta(:,end),'b')
hold off
title 'solution at right end as a function of time'; legend('alpha', 'beta');
grid on

function dy = fun(t,y,x)
  n = numel(y)/2;
  alpha = y(1:n);
  beta = y(n+1:2*n);
  dalpha = zeros(n,1);
  dbeta = zeros(n,1);
  dalpha(1) = alpha(1) + 2.2*beta(1);
  dalpha(2:n) = -0.6*(alpha(2:n)-alpha(1:n-1))./(x(2:n)-x(1:n-1));
  dbeta(1:n-1) = 1.32*(beta(2:n)-beta(1:n-1))./(x(2:n)-x(1:n-1));
  dbeta(n) = beta(n);
  dy = [dalpha;dbeta];
end

[Yidan]

Thanks! I will study this approach: method-of-lines approach.

Answer 2

[Bill Greene]

The reason for the error you are getting is that your boundary condition definition (p=0 and q=0) is invalid.

I modified your pdefun into an equivalent form and provided a valid boundary condition definition as shown below.

matlabAnswers_8_19_2024()

function matlabAnswers_8_19_2024
x = linspace(0,1,100);
t = linspace(0,1,30);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
if 0
    surf(x,t,u1)
    title('u_(x,t)')
    xlabel('Distance x')
    ylabel('Time t')
end

figure; plot(x, u1(end,:), x, u2(end,:));
title 'solution at final time'; legend('u1', 'u2');

figure; plot(t, u1(:,1), t, u2(:,1));
title 'solution at left end as a function of time'; 
legend('u1', 'u2');

figure; plot(t, u1(:,end), t, u2(:,end));
title 'solution at right end as a function of time'; 
legend('u1', 'u2');

end

function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve
c = [1; 1];
if 0
    f = [-0.6; 1.32].*u;
    s = [0; 0];
else
    f=[0 0]';
    s = [-0.6; 1.32].*dudx;
end
end

function u0 = pdeic(x) % Initial Conditions
u0 = [sin(pi*x); sin(pi*x)];
end

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
pl = [2.2*ul(2)+ ul(1); 0];
ql = [0; 1];
pr = [0; ur(2)];
qr = [1; 0];
end

[Yidan]

Sooooo cool!

Could you please tell me how you can modify the pdefun into an equivalent form? What dose "if 0" mean?

Could you please tell me if matrix p has elements equal to 0 in the boundary condition, does the element in the corresponding position q have to be 1 or some other value? What are the principles of boundary conditions?

[Bill Greene]

I simply moved the term you defined as f to the s term.

The "if 0" is just a trick for commenting out your original lines. This is just basic matlab code syntax that you should be familiar with.

[Yidan]

Thank you for your answer.

Could you please tell me if matrix p has elements equal to 0 in the boundary condition, does the element in the corresponding position q have to be 1 or some other value? What are the principles of boundary conditions?