Error in running the code for bvp4c with second order equation, how can I solve this issue?

Question

Hussein Al-Tamreh 2021-3-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/789539-error-in-running-the-code-for-bvp4c-with-second-order-equation-how-can-i-solve-this-issue

评论： Star Strider 2021-4-1

采纳的回答： Star Strider

在 MATLAB Online 中打开

I will try to explain briefly what happened, since i can't identify the location of my error.

I typed in my 2nd order equation in the "bvpode5" script as follows;

function dydx = bvpode5(x,y)
global Pe Theta K A ;
%
dydx = [y(2) 
  Pe(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
end

while the original equation is like this 1/Pe (d^2 α)/(dx^2 ) - dα/dx + θ[( α)/(K(A-α)+α)][A-α] = 0

my y(1) is α , while y(2) is dα/dx

As for my boundary condition's script is as follows:

function res = bvpbc5(ya,yb)
global Pe
res = [ya(1)-1-((1/Pe)*ya(2))
    yb(2)
    ];
end 

ya(1) represents α = 1 + (1/Pe)(dα/dx) at x = 0

ya(2) represents dα/dx = 0 at x = 1

Finally, the main script is as follows:

global Pe   Theta K  A;
% 
 % Kc = 0.482;  %   gCOD/gVSS
 Theta = 40;  %  day/day
 % Y = 0.212;   %   gVSS/gCOD
 Pe = 1000; 
 K = 0.123; 
 % Xo = 0.1;  % g/L
 % So = 0.1;  % g/L 
 A= 5.72;   % A = ((Xo/(Y*So))+1)
 xlow = 0;
 xhigh = 1; 
 y=[];
 solinit = bvpinit(linspace(xlow,xhigh,100),[1 0]);
 sol = bvp4c(@bvpode5,@bvpbc5,solinit);
 plot(sol.x,sol.y(1,:))
 xlabel('Reactor Length Z/L')
 ylabel('Substrate Fraction')
 axis ([0 1 0.2 1]);
%
fh=figure(1);
set(fh,'color','white')
title('Substrate Profile Across the Reactor')
disp (sol.y(1,end));

I tried to run the code from this script. However, it shows me this error, how can I solve this issue?;

MainpADM
Array indices must be positive integers or logical values.
Error in bvpode5 (line 6)
  Pe(y(2)-Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
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 MainpADM (line 16)
 sol = bvp4c(@bvpode5,@bvpbc5,solinit);
 

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

请先登录，再进行评论。

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

Answer 1

Star Strider 2021-4-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/789539-error-in-running-the-code-for-bvp4c-with-second-order-equation-how-can-i-solve-this-issue#answer_664014

在 MATLAB Online 中打开

There is a missing operator (most likely multiplication):

Pe(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
  ↑ ← HERE

This will probably work:

Pe*(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];

That should at least solve the problem you posted. I did not see any other problems with your code, however I did not run it.

It is always best to not use global variables. Pass the other parameters as additional parameters instead. See Passing Extra Parameters for an extended discussion.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2021-4-1

在 MATLAB Online 中打开

I am not certain what the problem is, however ‘sol.y(2,:)’ appears to be closer to the solution you want, instead of ‘sol.y(1,:)’.

Other than that, check the constants (parameters) and be certain the expreession is coded correctly, since it otherwise appears to do what you want.

The (slightly revised) code —

function dydx = bvpode5(x,y,Pe,Theta,K,A)
%
dydx = [y(2) 
        Pe*(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
end
function res = bvpbc5(ya,yb,Pe)
res = [ya(1)-1-((1/Pe)*ya(2))
       yb(2)];
end 
guess = @(x) [exp(x); -exp(-x)];
 % Kc = 0.482;  %   gCOD/gVSS
 Theta = 40;  %  day/day
 % Y = 0.212;   %   gVSS/gCOD
 Pe = 1000; 
 K = 0.123; 
 % Xo = 0.1;  % g/L
 % So = 0.1;  % g/L 
 A= 5.72;   % A = ((Xo/(Y*So))+1)
 xlow = 0;
 xhigh = 1; 
%  y=[];
 solinit = bvpinit(linspace(xlow,xhigh,150),guess);
 sol = bvp4c(@(x,y)bvpode5(x,y,Pe,Theta,K,A), @(ya,yb)bvpbc5(ya,yb,Pe),solinit);
  
 figure
 plot(sol.x,sol.y(2,:))
 xlabel('Reactor Length Z/L')
 ylabel('Substrate Fraction')
%  axis ([0 1 0.2 1]);
axis([0  1    -1  10])
%
fh=figure(1);
% set(fh,'color','white')
title('Substrate Profile Across the Reactor')
disp (sol.y(1,end));