pdepe - changing I.C 'midstream'

1 次查看（过去 30 天）
Eric
Eric 2015-7-29
编辑： Hashim 2021-4-15
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.

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2015-7-30
编辑：Torsten 2015-7-30
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.
  2 个评论
Eric
Eric 2015-7-30
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
Hashim 2021-4-14
编辑：Hashim 2021-4-15
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
%%
%%

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 PDE Solvers 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by