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:
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Torsten
2024-8-15
编辑:Torsten
2024-8-15
"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 个评论
Torsten
2024-8-19
编辑:Torsten
2024-8-20
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
更多回答(1 个)
Bill Greene
2024-8-19
编辑:Torsten
2024-8-19
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
3 个评论
Bill Greene
2024-8-19
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 PDE Solvers 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!