How to code for Relation boundary in pdepe
Afficher commentaires plus anciens
I have wrote the code using pdepe, but showing error, i cannot code for the realtion boundary condition and also it fails read parameters in boundary conditions. Please, correct the mistakes in the code.
function Daniel
m = 0;
x = linspace(0,1);
t = linspace(0,10000);
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2= sol(:,:,2);
u3= sol(:,:,3);
%figure
plot(x, u2(end,:),'.',...
'MarkerEdgeColor',[0.39 0.83 0.07],...
'LineWidth',2,'MarkerSize',8)
% plot(x,u1(end,:),'MarkerFaceColor',[1 0.6 0.6],'MarkerEdgeColor',[1 0 0],...
% 'MarkerSize',6,...
% 'Marker','diamond',...
% 'LineWidth',1,...
% 'LineStyle','none')
%legend('numerical')
%surf(x,t,u2)
hold on
%title('u1(x, t)')
%xlabel('Distance x')
%ylabel('u1(x,t)')
%figure
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(~,~,u,DuDx)
%alpha=1; beta=1; gamma=0.5; sigma=0.5; varphi=0.6; omega=3;mu=2;
alpha = 0.5;beta = 0.6;gamma = 0.7;sigma = 0.5;phi=0.8; mu=5; v=0.5;w=10;
c = [1;1;1];
f = [1;1;1].* DuDx;
F1=-phi^2.*u(1).*u(2)./(u(1).*u(2)+v.*u(2)+w.*u(1));
F2=-phi^2.*u(1).*u(2)./(u(1).*u(2)+v.*u(2)+w.*u(1));
F3=phi^2.*u(1).*u(2)./(u(1).*u(2)+v.*u(2)+w.*u(1));
s = [F1;F2;F3];
% --------------------------------------------------------------
function u0 = pdex4ic(~)
u0 = [1;1;1];
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(~,ul,~,ur,~)
pl = [0;ul(2);0];
ql = [1;0;0];
pr = [ur(1);ur(2);ur(3)];
qr = [-mu.*(1-l./alpha);sigma.*(1-m./be);gamma.*n];
Réponses (1)
You didn't supply the parameter values in pdex4bc for mu,l,alpha,sigma,m,be,gamma and n.
Further setting pl(3)=ql(3)=0 is not allowed. Setting C(0,t) = 0 gives pl(3)=ul(3),ql(3)=0. The boundary condition for B at 0 in your screenshot is strange and I don't understand it.
Your relations at x=1 are not coded correctly:
qr(1) = 1;
pr(1) = -mu1*(1-ur(1)/alpha1)
qr(2) = 1;
pr(2) = -mu2*(1-ur(2)/alpha2)
qr(3) = 1;
pr(3) = -mu3*ur(3);
Catégories
En savoir plus sur Programming dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!