How to use pdepe with arbitrary initial condition

I want to use pdepe with an arbitrary initial condition, rather than a defined symbolic equation.

The problem is that pdepe does not allow this, due to some stuff around lines 229-239 in pdepe.m where it uses single value from the x-vector to calculate a single value of the initial condition rather than just taking the assigned initial condition wholesale.

The code I want would look something like this

function sol = actual_pde(init_c)
%  Solving u_t  = d((u_x)^2)/dx 
%               = 2*u_xx * u_x ;
%
% FYI: init_c = sin(x)
x = linspace(0,1,20);
t = linspace(0,.1,5);
sol = pdepe(0,@pd_pde,@pd_initc,@pd_bc,x,t);

% --------------------------------------------------------------
function [c,f,s] = pd_pde(x,t,u,DuDx)
c = 1;
f = abs(DuDx)*DuDx;
s = 0;

% --------------------------------------------------------------
function u0 = pd_initc(init_c)
 u0 = init_c

% --------------------------------------------------------------
function [pl,ql,pr,qr] = pd_bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

but matlab will only accept

function u0 = pd_initc(x)
u0=sin(x);

It seems like the only solutions would be to edit the pdepe.m file around these lines:

% Form the initial values with a column of unknowns at each
% mesh point and then reshape into a column vector for dae.
temp = feval(pd_initc,xmesh(1),varargin{:});
if size(temp,1) ~= length(temp)
  error(message('MATLAB:pdepe:InvalidOutputICFUN'))
end

npde = length(temp);
y0 = zeros(npde,nx);
y0(:,1) = temp;
for j = 2:nx
  y0(:,j) = feval(pd_initc,xmesh(j),varargin{:});
end

But that seems like a huge workaround. Am I missing something obvious here?

Answer 1

[Jonathan]

Well, here is a solution. It involves using global variables to get around matlab's implementation of the pdepe algorithm.

Obviously this is not best practice but, hey, it works*

function sol = pde_solve()
%  Solving u_t  = d((u_x)^2)/dx 
%               = 2*u_xx * u_x ;

%Setup vars
x = linspace(0,1,20);
t = linspace(0,.1,5);

%define global vars (Ugh) to get around limitations in defining initial
%conditions
global x_vector;
global init_c_vector
init_c_vector = sin(x);
x_vector = x;

% Do solve
sol = pdepe(0,@pd_pde,@pd_initc,@pd_bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);

%%%%%%%%%%%%%%%%%%%
% Snipped out boring code
%%%%%%%%%%%%%%%%%%%

function u0 = pd_initc(x)
  %initialise global variables    
  global init_c_vector;
  global x_vector;
  % find the index of the x value pdepe wants to call
  [R,C] = find(x_vector==x,1,'first');
  %use the index to return the value of interest
  u0 = init_c_vector(R,C);

Any code not stated here is identical to that in the question. It is trivial to turn this into a solve for input x,t vectors

*Matlab 2014b, I guess it will stop working at some point if pdepe.m gets upgraded.

Answer 2

[Ojotule Onoja]

For multicomponent equations and components, how do you define initial conditions?

[Hafish Mahdi]

use matrix