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일
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개)

태그

제품


릴리스

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by