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

조회 수: 2(최근 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개)

범주

Find more on Test Model Components in Help Center and File Exchange

태그

제품


릴리스

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by