pdepe with cross elements in the differentiation "c"

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?

[Walter Roberson]

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.

[Francesca Turco]

Thanks, it reads much better now.

Answer 1

[Bill Greene]

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

[Francesca Turco]

This looks amazing! And it may also open up more options for my code - thank you so much! I would definitly like to try this tonight, and recode my equations to fit your new pde solver. In the link you refer to, I see a GitHub repository - I would like to read up a bit if there are any more instructions or similar. If there is any document that I shoudl download to read etc let me know. I'll try to parse this tonight and get back to you!

While I'm here :) I am also looking for a way to add an integral term to the second equation, ie make a PID feedback controller, but of course there is no way to access the previous time steps in the pde. I looked at your new solver pde1dm before, because it can add a ODE to the system, which should be able to access the internal x values. But I concluded that may not work for me since I need the previous t values, and my system would still have all PDEs, I can't transform it into a version with a PDE+ODE. But maybe I need to work more on that.

[Bill Greene]

Installation is straightforward. The easiest way is to download a zip file from github, unzip it to some directory, and then add that directory to your matlab path.

There is a pdf version of the user manual in the documents directory that covers the essentials of pde1dm. Unfortunately it is not nearly as complete or clear as I would like. If, after reading the document, you have questions about things that are unclear or simply missing, please let me know and I'll try to improve it!

[Francesca Turco]

I will try the installation as you suggest! Found the pdf, it seems really complete, I'll see if I can adapt that to the case you show here, and then to mine - I'll report back :)

Thank you so much - I'd also like to get ahead and reference your work and help on this. I can use the reference on the page you mentioned: Bill Greene (2023). pde1dM ( https://github.com/wgreene310/pde1dm/releases/tag/v1.3 ), GitHub. Retrieved August 29, 2023. Is there anything else I can add to the paper (that should go in at the end of the year)? We do use acknowledgements, on top of references, but I am not sure what you prefer.

thanks again, more when I parse this out!

[Bill Greene]

Referencing the github page is fine.

[carlos Hernando]

Hi @Bill Greene. Thank you for your code. I do have a quick question, is it possible to use vectorize calculation in your coupledMassExample.

Thank you!

Answer 2

[Torsten]

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.

[Francesca Turco]

In the general description of systems of PDEs, c is expressed as a matrix, diagonal in that specifica case - so I was trying to understand if non-diagnonal matrices were possible too. This is in order to address systems that have a dUdt(i) that depends on dUdt(i+n).

If this is not possible with pdepe, then it answers my question.

[Torsten]

Why not adding +fPd*U(3)*(right-hand side of the dU1/dt expression) + fPd*U(1)*(right-hand side of the dU3/dt expression) to the right-hand side of the dU2/dt expression ? Then your matrix C is the 5x5 identity matrix.

But be careful in the definition of pl, ql, pr and qr for U(2) then !

[Francesca Turco]

Very interesting idea, trying it out today, I'll report back!

The only issue is that dUdt1 and dUdt3 contain an f element, which has a second derivative (d/dxdUdx1) and d/dx(dUdx3), then mulpilied by U3 and U1 respectively. These have to be put in f, but U1*dUdx3 in f gets derived too, so I get dUdx1*dUdx3+U1*(dUdx3)2 etc. I am trying to find a way to subtract those extra elements that don't belong...

[Francesca Turco]

Ok, this shoudl work, but it gives an error...

I added the d/dx(dUdx*U) terms to f, and subtracted the 2*dUdx(1)*dUdx(3) term in s, which should give the right UdUdx2 in total. But it gives me the error of wrong size vectors:

"Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 5: Error in burn_TnP_piD (line 8)

[sol] = pdepe(m,@(x,t,U,dUdx)pdefun(x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPd,fPrad,Snbi,Spump, ..."

The full code is attached, the coefficients are passed from burn_inputs and Lzinputs, etc.

All the f, s and c entries seem to have the right format...

function [x,t,T,ne,Prad,PrZ,Pbrem,Pline,Pin,Te,Lz,Lzb,Palpha]=burn_TnP_piD;

x=linspace(0,1,151); t=linspace(0,2,100); makeplots=0; 
burn_inputs; load Lzinputs; makeplots=0; Te=newTe; Lz=newLzKr; m=1; 

Vol=2*pi^2*1.75*(x*0.56).^2*1.85*(1-(1-8/3/pi)*0.59/1.75); 

[sol] = pdepe(m,@(x,t,U,dUdx)pdefun(x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,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);
T=sol(:,:,1); Pin=sol(:,:,2); ne=sol(:,:,3);

clear nutty; nutty=sol(:,:,4);   % b_eff
for i=1:size(nutty,2), bE(:,i)=gradient(nutty(:,i),t); end
for i=1:length(t), Lzb(i,:)=bE(i,:)/1e34./cz0psm; end 
PrZ=bE.*ne.^2; Pbrem=brem*0.00534*2*ne.^2.*sqrt(T); Pline=liner*0.01*2*ne.^2; Prad=PrZ+Pbrem+Pline; 

clear nutty; nutty=sol(:,:,5);   % Palpha
for i=1:size(nutty,2), Palpha(:,i)=gradient(nutty(:,i),t); end

disp(['  T0=' num2str(T(end,1))]); 
disp(['  Palpha0=' num2str(Palpha(end,1))])
if makeplots == 1, plotburn_delay, end
end

function [c,f,s] = pdefun (x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPd,fPrad,Snbi,Spump, ...
                           xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift)
intT0=pchip(xx,T0p,x); intPin0=pchip(xx,Pin0p,x); intn0=pchip(xx,n0p,x); 
intnep=pchip(xx,neped,x); intcz=pchip(xx,cz0psm,x); intV=pchip(xx,Vol,x);
intDT=pchip(xx,DT,x); intDN=pchip(xx,DN,x); intTout=pchip(xx,Tout,x);
intfPp=pchip(xx,fPp,x); 

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=U(1)^2*(0.78*U(3))^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(U(1)+tshift)-0.0089152*(U(1)+tshift).^2+0.000103742*(U(1)+tshift).^3);
berror=U(1)*U(3)-intT0*intn0; 

c=[1 1 1 1 1]';
f=[intDT*dUdx(1) -fPd*(U(3)*dUdx(1)+U(1)*dUdx(3)) intDN*dUdx(3) 0 0]'; 
s=[fTpin*U(2)-fTpr*Pr-intTout*U(1)+alfa*Pa, ...  % Te
   -intfPp*(berror)+fPrad*Pr -fPd*U(3)*(fTpin*U(2)-fTpr*Pr-intTout*U(1)+alfa*Pa) -fPd*U(1)*(Snbi*U(2)+intnep-Spump*U(3)) +2*fPd*dUdx(1)*dUdx(3), ...    % 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 
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)]';
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
pr = [ur(1)-0.01 0 ur(3)-0.06 0 0]';
qr = [1 1 1 1 1]';
end

[Torsten]

Include

size(c)
size(f)
size(s)

to check whether c,f and s are column vectors of length 5.

Note that your boundary condition for U(2) becomes a problem. By setting

f(2) = -fPd*(U(3)*dUdx(1)+U(1)*dUdx(3))

it reads now

-fPd*(U(3)*dUdx(1)+U(1)*dUdx(3)) = 0

which most probably is not what you want to impose.

[Francesca Turco]

It works!! The error was because... I had 2 spaces between sum terms in the s(2) expression... which told s to have extra elements *eyeroll goes here* Even the BC seems to be ok: I am now using pl(2)=0, ql(2)=0, pr(2)=ur(2)-0.02 and qr(2)=0, which seems to be keeping the Pin (U(2)) down to the low value I need at the outer edge.

I had to put a differentiable function in fPd instead of a scalar, because I have profiles with a "pedestal", ie a relatively sharp drop near the outer edge, and the derivative of that was giving edge spikes in Pin. The effect of these new terms should be constrained near the centre of the spatial domain anyway, so getting a coefficient that goes to zero smothly was necessary anyway.

I tested this script with fPd=0 against the one that doesn't have all the new terms, just as a sanity check, and it gives exactly the same results! fPd (the derivative gain, essentially) needs to be kept in check because it can prevent convergence if it gets too big, but this also is a typical behaviour for this physics topic, I believe.

I can't thank you enough. What is the best way to reference your help and work? If I can make this work generally, it is supposed to go into a Phys. Rev Letters and/or a Nucl. Fusion paper late rin the year, I would like to acknowledge your participation there.

[Torsten]

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.

[Francesca Turco]

Agreed, I will try to implement Bill's solver tomorrow. The two ODE without spatial derivative here are just to extract two quantities that otherwise are difficult to check, ie they are conponents of the first 3 equations that I wanted to independently evaluate. I can actually remove both of them, I believe. The behaviour that I see when I scan the fPd vector of coefficients seems consistent with the origin profiles, so it is another indication that this method should be correct, besides recovering the results of the original system when these terms are set to zero (but exist in the system). I'll keep checking...