How to plot parameters within bvp4c
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Partial Differential Equation Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!