pdepe with cross elements in the differentiation "c"
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
x=linspace(0,1,151); t=linspace(0,10,100);
[sol] = pdepe(m,@(x,t,U,dUdx)pdefun(x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPi,fPd,fPrad,Snbi,Spump, ...
xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift), ...
@(x)pdeic(x,Te,Lz,xx,cz0psm,T0p,n0p,Pin0p,alfa,tshift,Vol), ...
@pdebc, ...
x,t);
function [c,f,s] = pdefun (x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPi,fPd,fPrad,Snbi,Spump, ...
xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift)
intT0=1; intPin0=1; intn0=1;
intnep=1; intcz=1; intV=1;
intDT=1; intDN=1; intTout=1;
intfPp=1; intfPi=1; intfPd=0.01;
Pr=L*1e34*interp1(Te,Lz,U(1),'pchip','extrap')*U(3)^2*intcz+ brem*0.00534*2*U(3)^2*sqrt(U(1))+ liner*0.01*2*U(3)^2;
Pa=0.01*U(3);
berror=U(1)*U(3)-intT0*intn0; inter=trapz(berror,xx);
%c=[1 1 1 1 1]';
%f=[intDT*dUdx(1) 0 intDN*dUdx(3) 0 0]';
%s=[fTpin*U(2)-fTpr*Pr-intTout*U(1) + alfa*Pa, ... % Te
% -intfPp*(berror)-intfPi*inter+fPrad*Pr, ... % Pin
% Snbi*U(2)+intnep-Spump*U(3), ... % ne
% 1e34*interp1(Te,Lz,U(1),'pchip','extrap')*intcz, ... % coefficient x Prad(Lz)
% alfa*Pa]'; % Palpha
c=[ 1 0 0 0 0 % U1 =T
-fPd*U(3) 1 -fPd*U(1) 0 0 % U2 =Pin
0 0 1 0 0 % U3 =ne
0 0 0 1 0 % U4 =Lz
0 0 0 0 1]; % U5 =Palpha
f=[intDT*dUdx(1) 0 0 0 0
0 0 0 0 0
0 0 intDN*dUdx(3) 0 0
0 0 0 0 0
0 0 0 0 0];
s=[fTpin*U(2)-fTpr*Pr-intTout*U(1)+alfa*Pa 0 0 0 0
0 -intfPp*(berror)-intfPi*inter+fPrad*Pr 0 0 0
0 0 Snbi*U(2)+intnep-Spump*U(3) 0 0
0 0 0 1e34*interp1(Te,Lz,U(1),'pchip','extrap')*intcz 0
0 0 0 0 alfa*Pa];
end
function u0 = pdeic(x,Te,Lz,xx,cz0psm,T0p,n0p,Pin0p,alfa,tshift,Vol)
intT0=pchip(xx,T0p,x); intPin0=pchip(xx,Pin0p,x); intn0=pchip(xx,n0p,x);
intcz=pchip(xx,cz0psm,x); intLz0=interp1(Te,Lz,intT0,'pchip','extrap'); intV=pchip(xx,Vol,x);
%u0=[intT0, intPin0, intn0, 1e34*intLz0*intcz, alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)]';
u0=[intT0 0 0 0 0
0 intPin0 0 0 0
0 0 intn0 0 0
0 0 0 1e34*intLz0*intcz 0
0 0 0 0 alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)];
u0=[intT0, intPin0, intn0, 1e34*intLz0*intcz, alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)]';
end
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
%pl = [0 0 0 0 0]'; % Te, Pin, ne, coeff, Palpha
%ql = [1 1 1 1 1]';
%pr = [ur(1)-0.01 ur(2)-0.02 ur(3)-0.06 0 0]'; % Te, Pin, ne, coeff, Palpha
%qr = [1 1 1 1 1]';
pl = zeros(5,5); % Te, Pin, ne, coeff, Palpha
ql = ones(5,5);
pr = [ur(1)-0.01 0 0 0 0
0 ur(2)-0.02 0 0 0
0 0 ur(3)-0.06 0 0
0 0 0 0 0
0 0 0 0 0]; % Te, Pin, ne, coeff, Palpha
qr = ones(5,5);
end
Hello, I have used the pdepe to solve a system of 5 coupled PDEs, that works fine. I want to add two terms in the LHS dUdt that makes "c" be a combination of dUdt(1)+dUdt(2)+dUdt(3). Is this possible? Or does the LHS of pdepe "c" can only have non-zero diagonal values?
More basically: can pdepe accept matrices or just vector inputs in c, s, and f?
The code below is the main script that I made. It works fine with the c, f, s coefficients that are commented out. With the c, f, s expressed as matrices, it gives the error "Error using pdepe: Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 5."
I am not sure why this is. Am I expressing the matrices in the wrong way, or is this just not possible with pdepe?
2 commentaires
Walter Roberson
le 28 Août 2023
Please edit your post to format the code.
Delete the current code, then press the > button in the CODE toolbar over the edit window. That will create a code insertion window. Copy the code from your editor and paste it in at that point, and it will be stored the same as is in your editor.
Réponses (2)
Bill Greene
le 29 Août 2023
If you are willing to consider an alternative to pdepe, I have written a PDE solver that mostly supports the same input syntax as pdepe but has several enhancements including the capability for a non-diagonal c-matrix. It can be downloaded here.
Here is a very simple two-equation PDE system that shows how you can define a non-diagonal c-matrix where the terms also depend on the dependent variables in the PDE.
function coupledMassExample
L=1;
n=21;
n2 = ceil(n/2);
x = linspace(0,L,n);
tf=.1;
nt=20;
t = linspace(0,tf,nt);
pdeFunc = @(x,t,u,DuDx) pdeDefinition(x,t,u,DuDx);
icFunc = @(x) icDefinition(x, L);
bcFunc = @(xl,ul,xr,ur,t) dirichletBC(xl,ul,xr,ur,t);
m=0;
sol = pde1dm(m, pdeFunc,icFunc,bcFunc,x,t);
u=sol(:,:,1);
v=sol(:,:,2);
figure; plot(t, u(:, n2), t, v(:, n2)); grid on;
xlabel 'time'
ylabel 'u,v'
title 'solution at center';
legend('u1', 'u2', 'Location', 'west')
figure; plot(x, u(end, :), x, v(end, :)); grid on;
xlabel 'x'
ylabel 'u,v'
title('solution at final time');
legend('u1', 'u2', 'Location', 'south')
fprintf('Solution: Time=%g, u_center=%g, v_center=%g\n', ...
t(end), u(end,n2), v(end,n2));
end
function [c,f,s] = pdeDefinition(x,t,u,DuDx)
c=[3-.1*u(2) 1+.2*u(1)
1+u(1) 2-.3*u(2)];
f = DuDx;
s=[0 0]';
end
function u0 = icDefinition(x,L)
u0 = [sin(pi*x/L); 3*sin(pi*x/L)];
end
function [pl,ql,pr,qr] = dirichletBC(xl,ul,xr,ur,t)
pl = ul;
ql = [0 0]';
pr = ur;
qr = [0 0]';
end
4 commentaires
Torsten
le 28 Août 2023
Déplacé(e) : Torsten
le 28 Août 2023
More basically: can pdepe accept matrices or just vector inputs in c, s, and f?
Why not reading the documentation of "pdepe" ?
pl, ql, pr and qr must be vectors.
s must be a vector.
c must be a diagonal matrix with elements either zero or positive.
f must be a vector.
8 commentaires
Torsten
le 30 Août 2023
Modifié(e) : Torsten
le 30 Août 2023
I'd check the results with a second PDE solver (e.g. @Bill Greene 's code). Your model contains at least two ODEs without spatial derivatives (U(4) and U(5)), and the conditions pl(2)=0, ql(2)=0 practically mean that no boundary condition is set. Strictly speaking, "pdepe" is not designed for this kind of system of differential equations. Despite your enthusiasm, I'd classify the results as at least necessary to be evaluated.
Voir également
Catégories
En savoir plus sur Eigenvalue Problems dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!