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:

Question

Yidan 2024-8-15，9:02

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2145439-i-am-trying-to-solve-the-system-of-coupled-partial-differential-equations-described-in-the-attachmen

评论： Yidan 2024-8-21，12:37

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')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2024-8-15，9:23

编辑：Torsten 2024-8-15，10:20

"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 .

6 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2024-8-19，13:23

编辑：Torsten 2024-8-20，0:23

在 MATLAB Online 中打开

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 2024-8-20，3:20

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

请先登录，再进行评论。

Answer 2

Bill Greene 2024-8-19，11:11

编辑：Torsten 2024-8-19，13:37

在 MATLAB Online 中打开

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

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Bill Greene 2024-8-19，12:24

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 2024-8-21，12:37

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?

请先登录，再进行评论。