MATLAB steady state (bvp4c) and transient (pdepe) plots do not match

The name of the pictureThe name of the pictureThe name of the pictureClash Royale CLAN TAG#URR8PPP











up vote
0
down vote

favorite












I'm trying to use MATLAB solve a reaction-diffusion equation
in both steady state (using bvp4c) and time dependency (using pdepe). When plotting both of the solutions, I do not get the same answer.
For my boundary conditions I have: at x = 0, C = 0 and at x = L, dC/dx = 0.



For my steady state case (bvp4c):



 function steady_state_bvp4c
 close all; clear all; clc;

 %Diffusion coefficient
 D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

 %Initial concentration at x=0
 L0 = 50;

 %Total length
 x_f = 20;

 %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 Receptors
 N_T =1.7;

 solinit = bvpinit(linspace(0,x_f,11),[0.5 0]);

 sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

 figure(1)
 plot(sol.x,(sol.y(1,:)),'LineWidth',1)
 title('Steady State')
 xlabel('Distance from Device (mum)')
 ylabel('Concentration (nM)')
 axis([0 x_f 0 L0])

 function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
 C_L = y(1);
 dC_Ldx = y(2);

 N_r = -(K_3.*K_4.*K_i.*(K_2 + C_L - ((8.*K_2.*C_L.^2.*N_T + 8.*K_2.*K_3.*C_L.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*C_L.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*C_L + 8.*K_2.*K_3.*K_i.*C_L.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*C_L.^2 + 4.*K_3.*C_L + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*C_L);
 N_R = (C_L./K_1).*N_r;
 N_r2 = (1./K_4).*(N_r).^2;
 N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
 N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
 N_R2 = ((C_L.^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.*C_L.*N_r))+...
 ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(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.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

 dy = zeros(2,1);
 dy(1) = dC_Ldx;
 dy(2) = (-R_L/D_ij);
 end

 function res = twobc(ya,yb)
 res = [ ya(1)-(L0); yb(2) ];
 end
 end


steady state plot



For my time dependent case:



 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 = 50; %c0 [nM]
 x_f =20; %Length of domain [um]
 maxt = 1; %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
 title('Time Dependent')
 xlabel('Distance from Device (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;
 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));

 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)
 % Boundary conditions for x = 0 and x = L;

 c0 = L0;

 % BCs: No flux boundary at the right boundary and constant concentration on
 % the left boundary

 % At x = 0, c0 = L0
 pl = ul-c0; 
 ql = 0; 

 % At x = L, dc/dx = 0
 pr = 0;
 qr = 1;

 end

 end


time dependent plot



Is there a difference in the way I am setting up both situations?
I would greatly appreciate any help!
Thank you!




differential-equations pde matlab




share|cite|improve this question



















  • There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway.
    – Jeff
    Jul 17 at 17:43















up vote
0
down vote

favorite












I'm trying to use MATLAB solve a reaction-diffusion equation
in both steady state (using bvp4c) and time dependency (using pdepe). When plotting both of the solutions, I do not get the same answer.
For my boundary conditions I have: at x = 0, C = 0 and at x = L, dC/dx = 0.



For my steady state case (bvp4c):



 function steady_state_bvp4c
 close all; clear all; clc;

 %Diffusion coefficient
 D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

 %Initial concentration at x=0
 L0 = 50;

 %Total length
 x_f = 20;

 %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 Receptors
 N_T =1.7;

 solinit = bvpinit(linspace(0,x_f,11),[0.5 0]);

 sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

 figure(1)
 plot(sol.x,(sol.y(1,:)),'LineWidth',1)
 title('Steady State')
 xlabel('Distance from Device (mum)')
 ylabel('Concentration (nM)')
 axis([0 x_f 0 L0])

 function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
 C_L = y(1);
 dC_Ldx = y(2);

 N_r = -(K_3.*K_4.*K_i.*(K_2 + C_L - ((8.*K_2.*C_L.^2.*N_T + 8.*K_2.*K_3.*C_L.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*C_L.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*C_L + 8.*K_2.*K_3.*K_i.*C_L.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*C_L.^2 + 4.*K_3.*C_L + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*C_L);
 N_R = (C_L./K_1).*N_r;
 N_r2 = (1./K_4).*(N_r).^2;
 N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
 N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
 N_R2 = ((C_L.^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.*C_L.*N_r))+...
 ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(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.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

 dy = zeros(2,1);
 dy(1) = dC_Ldx;
 dy(2) = (-R_L/D_ij);
 end

 function res = twobc(ya,yb)
 res = [ ya(1)-(L0); yb(2) ];
 end
 end


steady state plot



For my time dependent case:



 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 = 50; %c0 [nM]
 x_f =20; %Length of domain [um]
 maxt = 1; %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
 title('Time Dependent')
 xlabel('Distance from Device (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;
 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));

 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)
 % Boundary conditions for x = 0 and x = L;

 c0 = L0;

 % BCs: No flux boundary at the right boundary and constant concentration on
 % the left boundary

 % At x = 0, c0 = L0
 pl = ul-c0; 
 ql = 0; 

 % At x = L, dc/dx = 0
 pr = 0;
 qr = 1;

 end

 end


time dependent plot



Is there a difference in the way I am setting up both situations?
I would greatly appreciate any help!
Thank you!




differential-equations pde matlab




share|cite|improve this question



















  • There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway.
    – Jeff
    Jul 17 at 17:43













up vote
0
down vote

favorite









up vote
0
down vote

favorite











I'm trying to use MATLAB solve a reaction-diffusion equation
in both steady state (using bvp4c) and time dependency (using pdepe). When plotting both of the solutions, I do not get the same answer.
For my boundary conditions I have: at x = 0, C = 0 and at x = L, dC/dx = 0.



For my steady state case (bvp4c):



 function steady_state_bvp4c
 close all; clear all; clc;

 %Diffusion coefficient
 D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

 %Initial concentration at x=0
 L0 = 50;

 %Total length
 x_f = 20;

 %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 Receptors
 N_T =1.7;

 solinit = bvpinit(linspace(0,x_f,11),[0.5 0]);

 sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

 figure(1)
 plot(sol.x,(sol.y(1,:)),'LineWidth',1)
 title('Steady State')
 xlabel('Distance from Device (mum)')
 ylabel('Concentration (nM)')
 axis([0 x_f 0 L0])

 function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
 C_L = y(1);
 dC_Ldx = y(2);

 N_r = -(K_3.*K_4.*K_i.*(K_2 + C_L - ((8.*K_2.*C_L.^2.*N_T + 8.*K_2.*K_3.*C_L.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*C_L.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*C_L + 8.*K_2.*K_3.*K_i.*C_L.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*C_L.^2 + 4.*K_3.*C_L + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*C_L);
 N_R = (C_L./K_1).*N_r;
 N_r2 = (1./K_4).*(N_r).^2;
 N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
 N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
 N_R2 = ((C_L.^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.*C_L.*N_r))+...
 ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(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.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

 dy = zeros(2,1);
 dy(1) = dC_Ldx;
 dy(2) = (-R_L/D_ij);
 end

 function res = twobc(ya,yb)
 res = [ ya(1)-(L0); yb(2) ];
 end
 end


steady state plot



For my time dependent case:



 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 = 50; %c0 [nM]
 x_f =20; %Length of domain [um]
 maxt = 1; %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
 title('Time Dependent')
 xlabel('Distance from Device (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;
 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));

 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)
 % Boundary conditions for x = 0 and x = L;

 c0 = L0;

 % BCs: No flux boundary at the right boundary and constant concentration on
 % the left boundary

 % At x = 0, c0 = L0
 pl = ul-c0; 
 ql = 0; 

 % At x = L, dc/dx = 0
 pr = 0;
 qr = 1;

 end

 end


time dependent plot



Is there a difference in the way I am setting up both situations?
I would greatly appreciate any help!
Thank you!




differential-equations pde matlab




share|cite|improve this question











I'm trying to use MATLAB solve a reaction-diffusion equation
in both steady state (using bvp4c) and time dependency (using pdepe). When plotting both of the solutions, I do not get the same answer.
For my boundary conditions I have: at x = 0, C = 0 and at x = L, dC/dx = 0.



For my steady state case (bvp4c):



 function steady_state_bvp4c
 close all; clear all; clc;

 %Diffusion coefficient
 D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

 %Initial concentration at x=0
 L0 = 50;

 %Total length
 x_f = 20;

 %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 Receptors
 N_T =1.7;

 solinit = bvpinit(linspace(0,x_f,11),[0.5 0]);

 sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

 figure(1)
 plot(sol.x,(sol.y(1,:)),'LineWidth',1)
 title('Steady State')
 xlabel('Distance from Device (mum)')
 ylabel('Concentration (nM)')
 axis([0 x_f 0 L0])

 function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
 C_L = y(1);
 dC_Ldx = y(2);

 N_r = -(K_3.*K_4.*K_i.*(K_2 + C_L - ((8.*K_2.*C_L.^2.*N_T + 8.*K_2.*K_3.*C_L.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*C_L.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*C_L + 8.*K_2.*K_3.*K_i.*C_L.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*C_L.^2 + 4.*K_3.*C_L + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*C_L);
 N_R = (C_L./K_1).*N_r;
 N_r2 = (1./K_4).*(N_r).^2;
 N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
 N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
 N_R2 = ((C_L.^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.*C_L.*N_r))+...
 ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(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.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

 dy = zeros(2,1);
 dy(1) = dC_Ldx;
 dy(2) = (-R_L/D_ij);
 end

 function res = twobc(ya,yb)
 res = [ ya(1)-(L0); yb(2) ];
 end
 end


steady state plot



For my time dependent case:



 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 = 50; %c0 [nM]
 x_f =20; %Length of domain [um]
 maxt = 1; %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
 title('Time Dependent')
 xlabel('Distance from Device (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;
 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));

 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)
 % Boundary conditions for x = 0 and x = L;

 c0 = L0;

 % BCs: No flux boundary at the right boundary and constant concentration on
 % the left boundary

 % At x = 0, c0 = L0
 pl = ul-c0; 
 ql = 0; 

 % At x = L, dc/dx = 0
 pr = 0;
 qr = 1;

 end

 end


time dependent plot



Is there a difference in the way I am setting up both situations?
I would greatly appreciate any help!
Thank you!





differential-equations pde matlab





share|cite|improve this question










share|cite|improve this question




share|cite|improve this question









asked Jul 16 at 21:02









BaiSango

11




11











  • There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway.
    – Jeff
    Jul 17 at 17:43

















  • There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway.
    – Jeff
    Jul 17 at 17:43
















There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway.
– Jeff
Jul 17 at 17:43





There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway.
– Jeff
Jul 17 at 17:43
















active

oldest

votes











Your Answer




StackExchange.ifUsing("editor", function ()
return StackExchange.using("mathjaxEditing", function ()
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
);
);
, "mathjax-editing");

StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "69"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);

else
createEditor();

);

function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
convertImagesToLinks: true,
noModals: false,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);



);








 

draft saved


draft discarded


















StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2853853%2fmatlab-steady-state-bvp4c-and-transient-pdepe-plots-do-not-match%23new-answer', 'question_page');

);

Post as a guest






















active

oldest

votes













active

oldest

votes









active

oldest

votes






active

oldest

votes










 

draft saved


draft discarded


























 


draft saved


draft discarded














StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2853853%2fmatlab-steady-state-bvp4c-and-transient-pdepe-plots-do-not-match%23new-answer', 'question_page');

);

Post as a guest
































































Comments

Post a Comment

Popular posts from this blog

What is the equation of a 3D cone with generalised tilt?

Image
Clash Royale CLAN TAG #URR8PPP up vote 2 down vote favorite What is the equation of a 3D cone with generalized tilt? I've noticed that in most equations given to represent a cone, there is no parameter which defines the tilt of the cone in 3D space and that most of them have their point at the origin $(0,0,0)$ - I was wondering if anyone could give me a more generalized cone equation for a cone in any position in 3D space 3d share | cite | improve this question edited Jul 25 '16 at 0:05 ervx 9,425 3 13 36 asked Jul 24 '16 at 15:38 Charlene 11 2 2 Welcome to Math.SE! Could you please explain a few things to help people give an answer useful to you: 1. What do you mean by "generalized tilt"? 2. Are you asking about a right circular cone? Does the angle at the vertex matter? 3. What form of answer (implicit, parametric...?) are you looking for? 4. If this is homework, what tools do you have available, an...
Read more

Color the edges and diagonals of a regular polygon

Image
Clash Royale CLAN TAG #URR8PPP up vote 5 down vote favorite 1 Here is the problem: For what $n$ is it possible to color the edges and diagonals of an $n$-side regular polygon with $dfracbinomn23$ colors, such that you use every color exactly three times and for every color the three segments (edges or diagonals) with that color form a triangle? Trivially $nequiv 0 text or 1 pmod3$. I can also prove that if the statement is true for $k$ than it is true for $3k$ as well. How to finish? Please help! Thanks combinatorial-geometry share | cite | improve this question edited Aug 6 at 21:57 alcana 118 4 asked Aug 6 at 21:15 Leo Gardner 356 11 Even assuming "polyhpn" in the title is a typo for polygon , it is unclear what you mean by "the edges and sides". Aren't the side of a regular polygon the same as the edges? A regular $n$-sided polygon has $n$ edges. – hardmath Aug 6 at 21:20 3 $n$ ne...
Read more

Relationship between determinant of matrix and determinant of adjoint?

Image
Clash Royale CLAN TAG #URR8PPP up vote 2 down vote favorite We are studying adjoints in class, and I was curious if there is a relationship between the determinant of matrix A, and the determinant of the adjoint of matrix A? I assume there would be a relationship because finding the adjoint requires creating a cofactor matrix and then transposing it. linear-algebra determinant adjoint-operators share | cite | improve this question asked Nov 17 '17 at 17:32 CluelessCoder 32 5 For a square matrix $A$ of order $n$, we have $A(operatornameadj A)=(operatornameadj A)A=|A|I_n$ where $I_n$ is the identity matrix of order $n$. This should tell you about the determinant of the adjoint in terms of that of $A$ (use multiplicative property and other elementary properties of determinants). – Prasun Biswas Nov 17 '17 at 17:35 1 @CluelessCoder: see here: math.stackexchange.com/questions/516127/…...
Read more