bvp4c setup for second order ODE

Question

Jennifer Yang 2018-6-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406085-bvp4c-setup-for-second-order-ode

评论： 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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2018-6-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406085-bvp4c-setup-for-second-order-ode#answer_325062

在 MATLAB Online 中打开

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更早的评论隐藏 7更早的评论

Torsten 2018-6-19

在 MATLAB Online 中打开

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