pdepe - changing I.C 'midstream'

I'm using pdepe to solve a system of two nonlinear reaction-diffusion equations, Neumann B.C. At some time, t1, I would like to have pdepe use the solution, u(x,t1), a size(t)-by-size(x)-by-2 array, as the new Initial Conditions and have pdepe continue the solution from there. I want to do this because at t1 I want to change the values of some of the parameters in the reaction function specification (in pdefun) and get to the final solution. For the life of me, I cannot figure out how to pass u(x,t1) to icfun so that pdepe will use u(x,t1) as the I.C. I tried the anonymous function approach, but pdepe didn't seem to like that...perhaps I wasn't coding it correctly.

Answer 1

[Torsten]

global counter
counter=0;
...
solold=pdepe(m,@pdefunold,@icfunold,@bcfunold,xmesh,tspanold);
uold(1,:) = solold(end,:,1);
uold(2,:) = solold(end,:,2);
tspannew=...;
solnew=pdepe(m,@pdefunnew,@(x)icfunnew(x,uold),@bcfunnew,xmesh,tspannew);
...

function u0 = icfunnew(x,uold)
global counter
counter=counter+1;
u0(1)=uold(1,counter);
u0(2)=uold(2,counter);

Maybe you will have to adjust the variable "counter" from 0 to 1 in the calling program.

Best wishes

Torsten.

[Eric]

That solved it! The one additional thing I had to do was 'convert' the 2 u0 values to a column vector. Thank you sooooo much, and kudos for all the help you've given posters in the past!!!

[Hashim]

Hi, I think I am trying to do something similar where I have different IC for two different time spans. Would it be possible if I could get some help with my problem? @Torsten

function pdepe_philip_v3a

% 1-b Extend the above calculation to model a potential step to any chosen 
% potential assuming that the ratio of the surface concentrations of Os2+ 
% and Os3+ obeys the Nernst equation both before and after the potential
% step. Now we want to further add a potential step which will make it a
% double potential step. 

global D n F R Os_bulk 
n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D          = 5e-06;                              % cm^2/s 
Area       = 1;                                  % cm^2 
Os_bulk    = 1e-05;                              % mol/cm^3

N          = 10; 

m          = 0;                                  % Cartesian Co-ordinates

% logarithmized xmesh for semi infinite boundary condition 
x          = logspace(log10(0.000006), log10(0.06), N); x(1) = 0;


% Time Span-1
t1         = linspace(0, 10, N);                 % s      
c0         = Os_bulk; 
E          = 0.8; 
ic_arg     = @(x)c0.*ones(size(x)); 
sol1       = pdepe(m, @pdev3apde, @(x)pdev3aic(x,ic_arg),@(xl, ul, xr, ur, t)...
                   pdev3abc(xl, ul, xr, ur, t, E, c0), x, t1);
c1         = sol1(:, :, 1); 
I_num      = ( n*F*Area) * c1(1,:) .* sqrt(D./(t1*pi));  

% I Flux For Time Span-1
figure(1); 
plot(t1, I_num, 'r--', 'LineWidth',1);
xlabel('t (s)'); ylabel('I (A/cm^2s)'); 
hold on 

% Time Span-2
t2         = linspace(10, 20, N);                % s       
E          = 1.6; 
c0         = (c1(1,:))'; 
ic_arg     = @(x)(c0).*ones(size(x));
sol2       = pdepe(m, @pdev3apde, @(x)pdev3aic(x,ic_arg),@(xl, ul, xr, ur, t)...
                   pdev3abc(xl, ul, xr, ur, t, E, c0), x, t2);
c2         = sol2(:, :, 1); 
I_num1     = ( n*F*Area) * c2(1,:) .* sqrt(D./(t2*pi));  

% I Flux For Time Span-2
figure(1); 
plot(t2, I_num1, 'b--', 'LineWidth',1);
xlabel('t (s)'); ylabel('I (A/cm^2s)'); 
hold on 

end

%% pdepe Function That Calls Children
%%

function [a, f, s] = pdev3apde(x, t, u, DuDx)

global D

a           = 1; 
f           = -D*DuDx;
s           = 0;

end

%% Initial Condition
%%

function u0 = pdev3aic(x, ic_arg)

u0          = ic_arg(x);

end

%% Boundry Conditions 
%%

function [pl, ql, pr, qr] = pdev3abc(xl, ul, xr, ur, t, E, c0)

global n F R Os_bulk 
 
E0          = 0.2;
alpha       = exp((E-E0).*((n*F)/(R*298.15))); 
pl          = ul - (c0./(1+alpha));
ql          = 0;
pr          = ur - Os_bulk;
qr          = 0;     

end

%%
%%