Empirical source terms with PDEPE

Suppose I'm solving an equation:

With initial condition and boundary conditions and

Now I want to import the source term, as data, so I will have a Nx2 array [t Q] with data. Can I just treat that as I would a constant?

The code I have is:

function sol=temp_pde(t,R,X,source)
m=1; %Sets the geometry to cylindrical
global theta kappa h Q;
theta=X(1); kappa=X(2); h=X(3);
r=linspace(0,R,800);
t=source(:,1);
Q=source(:,2);
sol=pdepe(m,pdefun,icfun,bcfun,r,t);
end

function [c,f,s] = pdefun(r,t,u,dudx)
global Q theta kappa;
c = theta;
f = kappa*dudx;
s = Q;
end

function u0 = icfun(r)
u0 = 0;
end

function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
global h kappa;
pl = 0;
ql = 1;
pr = h*ur;
qr = kappa;
end

Where I obtain my t parameter from my data in source. Would this work?

[Bjorn Gustavsson]

Is a time-variable heating-rate that is the same over the entire region?

It is better to wrap the heating-rate into a function to make it independent of time, something like this:

Qfcn = @(t,r,u) ones(size(r))*interp1(source(:,1),source(:,2),'pchip');

Now will be a function of time that pdepe can evaluate at any time - it will use some time-steps (adaptive or fixed) that will not necessarily line up with the time-steps of your source.

(Check the documentation for how the right-hand-side function should be defined in terms of input variables and the like.)

HTH

[Matthew Hunt]

\dot{Q} is purely a function of time, so is homogeneous.

Answer 1

[Bjorn Gustavsson]

To the best of my understanding you should be able to define your pde-function something like this:

function [c,f,s] = pde_interpS(x,t,u,dudx,tQ,Q,theta,kappa)
c = theta;
f = kappa*dudx;
s = interp1(tQ,Q,t,'pchip');
end

Then you can convert that function to a function of t, x, u and dudx with the standard @()-trick. My preference is to avoid globals as much as possible, but you can use that too.

HTH

[Matthew Hunt]

Why would I need to use interp when I'm getting the time spacing from the data?

[Bjorn Gustavsson]

You get that time-spacing - but pdepe doesn't necessarily take only those timesteps. Just as the ODE-integrating functions are adaptive - i.e. take variable time-steps to reach the required accuracy. If you check the documentation you'll see:

The time integration is done with the ode15s solver. pdepe exploits the capabilities of ode15s for solving the differential-algebraic equations that arise when the PDE contains elliptic equations, and for handling Jacobians with a specified sparsity pattern.

HTH

[Matthew Hunt]

However, this isn't where I get errors oddly. I get a problem with this line:

f = kappa*dudx;

I have absolutely no idea why this is causing errors. Any idea?

[Bjorn Gustavsson]

What's the error-message?

Have you tried to run the function with dbstop set to if error?

dbstop if error

That's the best way to figure out what's going on. I got something similar to work without too much of a problem.

HTH

[Matthew Hunt]

I spotted my error, I should have "@" signs, now I have a new problem, I get the following error:

246 [c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

247 if any([size(c,1),size(f,1),size(s,1)]~=npde)

248 error(message('MATLAB:pdepe:UnexpectedOutputPDEFUN',sprintf('%d',npde)))

Error using pdepe (line 248)

Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 1.

[Bjorn Gustavsson]

Seems like your pde-function that should return your c, f and s either doesn't return all three of those variables, or some of them are not set to column vectors - either one or more of them are empty or of the wrong size - perhaps one of them becomes a row-vector. If you set a debug-stop in the pde-function you can check the sizes of the outputs there.

HTH