How to plot parameters within bvp4c
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (0)
Voir également
Catégories
En savoir plus sur Partial Differential Equation Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!