Index exceeds array bounds error message in a three-point bvp

조회 수: 1 (최근 30일)
Dursman Mchabe
Dursman Mchabe 2018년 8월 28일
Hi everyone, I'm trying to solve a three-point bvp, following:
https://www.mathworks.com/help/releases/R2018a/matlab/ref/bvpinit.html
However, I get the error message:
"Index exceeds array bounds.
Error in ConcentrationProfilesinLimestoneSlurryLiquidFilm>materialbalanceodes (line 63) dydx(1) = y(6);
Error in bvparguments (line 105) testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 130) bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in ConcentrationProfilesinLimestoneSlurryLiquidFilm (line 4) sol = bvp4c(@materialbalanceodes,@materialbalancebc,solinit);"
What can I try? Please see the code on the attachment or below:
x = [0, 10e-6, 20e-6, 20e-6, 30e-6, 40e-6];
solinit = bvpinit(x,[4.03e-08, 3.39e-02, 2.60e-01, 5.45e-01, 4.82e-03, 0, 0, 0]); sol = bvp4c(@materialbalanceodes,@materialbalancebc,solinit);
y = deval(sol,x); plot(x,y(1,:),x,y(2,:),x,y(3,:),x,y(4,:),x,y(5,:)) legend('SO_2','HSO_{3}^{-}','SO_{3}^{2-}','HCO_{3}^{-}','CO_{3}^{2-}')
global A B C D E F G H I J K L M N O P Q R S T U
%PARAMETERS A = sym(9.4e-4); % ks = solid-liquid mass transfer (m/s) B = sym(3.69e2); % Ap = solid-liquid surface area (m2/m3) C = sym(2.89e-9); % DSO2 = diffusivity of SO2 (m2/s) D = sym(1.6e-9); % DCO32 = diffusivity of CO32- (m2/s) E = sym(4.92e-5); % CBs = concentration of Ca2+ on limestone surface (mol/m3) F = sym(1.6e-8); % DH = diffusivity of H+ (m2/s) G = sym(6.24); % KSO2 = SO2 dissociation equilibrium (mol/m3) H = sym(2.35e-9); % DHSO3 = diffusivity of HSO3 (m2/s) I = sym(5.68e-5); % KHSO3 = HSO3 dissociation equilibrium constant (mol/m3) J = sym(1.69e-9); % DSO32 = diffusivity of SO32 (m2/s) K = sym(149); % HSO2 = Henry's constant (m3Pa/mol) L = sym(47.28325); % PSO2 = partial pressure of SO2 (Pa) M = sym(4.14e-6); % kga = gas-side mass transfer product (1/s) N = sym(8.314); % R = universal gas constant (J/mol.K) O = sym(323.15); % T = reactor tempeature (K) P = sym(2.09); % CH = concentration of H+ at the gas-liquid interface (mol/m3) Q = sym(9.333E-06); % CH = concentration of H+ at the in the bulk slurry (mol/m3) R = sym(0.679); % S(IV) = total sulfur concentration in the bulk slurry (mol/m3) S = sym(4.88e-6); % C(IV) = total sulfur concentration in the bulk slurry (mol/m3) T = sym(1.7e-3); % KCO2 = CO2 dissociation equilibrium constant (mol/m3) U = sym(6.55e-8); % KHCO3 = HCO3 dissociation equilibrium constant (mol/m3)
% VARIABLES
% y(1) = CSO2 % y(2) = CHSO3 % y(3) = CSO32 % y(4) = CHCO32 % y(5) = CCO32
% EQUATIONS
% Region I
% dy(1)/dx = y(6); % d^2y(1)/dx^2 = (2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) / D*E)*E; % dy(2)/dx = y(7); % d^2y(2)/dx^2 = -(2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) /D*E)*E;
%Region II
% dy(3)/dx = y(8) % dy(4)/dx = y(4); % d^2y(3)/dx^2 = -(A*B / D)*(1 + (F*I*(y(3) / y(2)))/(D*E))*E; % dy(5)/dx = y(5)
function dydx = materialbalanceodes(y,x, A ,B, C, D, E, F, G, H, I)
dydx = zeros(8,1);
dydx(1) = y(6);
dydx(2) = (2 * A * B / C) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G )* (y(2)/y(1)) / (D*E)))*E;
dydx(3) = y(7);
dydx(4) = -((2 * A * B) / (H)) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G) * (y(2) / y(1)) / (D * E))) * E;
dydx(5) = y(8);
dydx(6) = y(4);
dydx(7) = -((A * B) / D) * (1 + (F * I * (y(3) / y(2)))/(D * E)) * E;
dydx(8) = y(5);
end
function res = materialbalancebc(ya,yb, J, K, L, M, N, O, P, Q) res = [(ya(1)- (K * L)); (ya(2)-((I * K * L) / P)); (ya(6)+ M * ((L /( N * O))- (K * L))); (yb(2) - ((R * G * Q)/(Q.^2 + G * Q + G * I))); (yb(3) - ((R * G * I)/(Q.^2 + G * Q + G * I))) ; yb(4) - ((S * T * Q)/(Q.^2 + T * Q + T * U))]; end

답변 (0개)

카테고리

Help CenterFile Exchange에서 Numeric Solvers에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by