bvp4c setup for second order ODE

18 次查看（过去 30 天）
Jennifer Yang
Jennifer Yang 2018-6-18
评论： Utkarsh Tiwari 2020-10-22
Hello, I'm trying to solve for this second order ODE in steady state using bvp4c with the boundary conditions where at x=0, C_L=1 and x=100, C_L=0.
function bvp4c_test
%Diffusion coefficient
D_ij= 1*10^-6;
%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,10,11),[1 0]);
sol = bvp4c(@odefcn,@twobc,solinit);
plot(sol.x,sol.y(1,:))
function dC_Ldx = odefcn(C_L,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)
dC_Ldx = zeros(2,1);
dC_Ldx(1) = C_L(2);
N_r = (-(1+(C_L./K_2))+...
sqrt(((1+(C_L./K_2)).^2)-4.*((2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3))))).*(-N_T)))./...
(2.*(2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3)))));
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));
dC_Ldx(2) = -R_L./D_ij;
end
function res = twobc(ya,yb)
res = [ ya(1)-1; yb(1) ];
end
end
When I run it, I encounter the error stating Attempted to access C_L(2); index out of bounds because numel(C_L)=1.
| | _Error in bvp4c_test/odefcn dC_Ldx(1) = C_L(2);
Error in bvparguments testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c [n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] = ...
Error in bvp4c_test sol = bvp4c(@odefcn,@twobc,solinit);_||
Where am I going wrong? I am still very new to MATLAB and appreciate any help that I can get. Thank you!

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2018-6-18
sol = bvp4c(@(x)odefcn(x,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);
  9 个评论
显示 7更早的评论
Torsten
Torsten 2018-6-19
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 = (-(1+(C_L./K_2))+...
sqrt(((1+(C_L./K_2)).^2)-4.*((2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3))))).*(-N_T)))./...
(2.*(2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3)))));
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
Utkarsh Tiwari
Utkarsh Tiwari 2020-10-22
Wow this was incredibly helpful as I faced a similar problem. Thank you!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by