How to pass solution of pdefun1 as IC-2 for pdefun2

1 次查看（过去 30 天）
Hashim
Hashim 2021-4-16
回答： Hashim 2021-4-17
my pde is
Given below is my effort so far...
function pdepe_philip_v3b
% 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
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
E = 0.8;
sol1 = pdepe(m, @pdev3bpde, @(x)pdev3bic(x,Os_bulk), @(xl, ul, xr, ur, t)...
pdev3bbc(xl, ul, xr, ur, t, E, Os_bulk), 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
c_new = (c1(N,:))';
E = -1.6;
sol2 = pdepe(m, @pdev3bpdeNEW, @(x)pdev3bicNEW(x,c_new),@(xl, ul, xr, ur, t)...
pdev3bbcNEW(xl, ul, xr, ur, t, E, c_new), x, t2);
end
%% pdepe Function
%%
function [a, f, s] = pdev3bpde(x, t, u, DuDx)
global D
a = 1;
f = D*DuDx;
s = 0;
end
%% pdepe Function
%%
function [a, f, s] = pdev3bpdeNEW(x, t2, u, DuDx)
global D
a = 1;
f = D*DuDx;
s = 0;
end
%% Initial Condition-1
%%
function u0 = pdev3bic(x, Os_bulk)
u0 = Os_bulk;
end
%% Initial Condition-2
%%
function u0 = pdev3bicNEW(x, c_new)
u0 = c_new.*size(x);
end
%% Boundry Condition-1
%%
function [pl, ql, pr, qr] = pdev3bbc(xl, ul, xr, ur, t, E, Os_bulk)
global n F R
E0 = 0.2;
alpha = exp((E-E0).*((n*F)/(R*298.15)));
pl = ul - (Os_bulk./(1+alpha));
ql = 0;
pr = ur - Os_bulk;
qr = 0;
end
%% Boundry Condition-2
%%
function [pl, ql, pr, qr] = pdev3bbcNEW(xl, ul, xr, ur, t, E, Os_bulk)
global n F R
E0 = 0.2;
alpha = exp((E-E0).*((n*F)/(R*298.15)));
pl = ul - (Os_bulk./(1+alpha));
ql = 0;
pr = ur - Os_bulk;
qr = 0;
end
%%
%%
  2 个评论
Hashim
Hashim 2021-4-16
编辑：Hashim 2021-4-16
It does not its just that previously i've been trying to pass the ic as a fucntion to icfun. In this attempt i've tried to make two instances of icfun one old and one "NEW". Otherwise its the same problem.
Also, I think this poster's problem is similar to mine in terms of icfun implementation but I can't seem to figure it out for the life of me.

请先登录，再进行评论。

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

采纳的回答

Hashim
Hashim 2021-4-17
Just so that anybody's thinking there is a way to pass solution from pdefun1 to icfunc2... here is a solution from @Torsten on this post. I am copying the code below (it is by Torsten of course).
function main
m = 0;
xmesh = linspace(0,20,101);
tspan = linspace(0,365*5*86400,201);
u0 = 280.0;
icarg = @(x) 0.01*ones(size(x));
sol = pdepe(m,@pdefun,@(x)icfun(x,icarg),@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,u0),xmesh,tspan);
w = sol(end,:);
plot(xmesh,w)
tspan2 = linspace(tspan(end),365*20*86400,201);
u0 = 0.0;
icarg = @(x)interp1(xmesh,w,x);
sol2 = pdepe(m,@pdefun,@(x)icfun(x,icarg),@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,u0),xmesh,tspan2);
w2 = sol2(1,:);
hold on
plot(xmesh,w2)
end
function [c,f,s] = pdefun(xmesh,tspan,u,DuDx)
c = 2;
f = 1e-8*DuDx;
s = -1e-7*u;
end
function u = icfun(x,icarg)
u = icarg(x);
end
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t,u0)
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Characters and Strings 的更多信息

标签

产品


版本

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by