How to plot parameters within bvp4c

1 次查看(过去 30 天)
Hello,
I'm solving for a reaction diffusion equation. I followed examples on how to plot the solution of differential equation and was able to do that. However, I'm not too sure how to plot the N_r, N_R, N_r2, N_rR, N_pR, N_R2, R_L relative to C_L or x.
function time_dependent_pdepe
clear all; close all; clc;
D_ij = 30; %Diffusion coefficient D (3.0*10^-7 cm^2/s -> 30 um^2/s)
L0 = 10; %c0 [nM]
x_f =40; %Length of domain [um]
maxt = 10; %Max simulation time [s]
m = 0; %Parameter corresponding to the symmetry of the problem
x = linspace(0,x_f,100); %xmesh
t = linspace(0,maxt,100); %tspan
sol = pdepe(m,@DiffusionPDEfun,@DiffusionICfun,@DiffusionBCfun,x,t,[]);
u = sol;
% Plotting
hold all
for n = linspace(1,length(t),10)
plot(x,sol(n,:),'LineWidth',2)
end
xlabel('Distance (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])
function [c,f,s] = DiffusionPDEfun(x,t,u,dudx)
D = D_ij;
%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;
%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;
%Total Recptors
N_T = 1.7;
N_r = -(K_3.*K_4.*K_i.*(K_2 + u - ((8.*K_2.*u.^2.*N_T + 8.*K_2.*K_3.*u.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*u.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*u + 8.*K_2.*K_3.*K_i.*u.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*u.^2 + 4.*K_3.*u + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*u);
N_R = (u./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (u./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (u./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((u.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;
R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*u.*N_r))+...
((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*u.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
(d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*u.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));
% PDE
c = 1;
f = D_ij.*dudx;
s = R_L;
end
function u0 = DiffusionICfun(x)
u0 = 0;
end
function [pl,ql,pr,qr] = DiffusionBCfun(xl,ul,xr,ur,t)
c0 = L0;
% At x = 0, c0 = L0
pl = ul-c0;
ql = 0;
% At x = L, dc/dx = 0
pr = 0;
qr = 1;
end
end
When I try plotting it, I get a blank plot. How would you go about plotting those values?
Thank you!

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by