Error in bvp4c code

Hello all!

I'm a research scholar and I'm trying to solve a coupled system of ODEs with boundary conditions. I'm using bvp4c solver. However, I'm getting the error, 'Index exceeds matrix dimensions' when I run the code. I'm attaching the equations that I'm solving as pdf. I'm also giving the code below. Please advice on how to proceed.

Thank You

R=1;
ya=0; yRc=0.86; yRp=0.91; yRb=1;
M=1.5;
alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
Rp=0.91;
Rb=0.95;
Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;
ud=P0*Da;

eps=1.0000e-04;
r=linspace(eps,R,100);
solinit = bvpinit(linspace(eps,R,100),[1,1,1,1,1,1,1,1,1]); 
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)base(r,y,P0,mu,hm,Rc,We,m,Da,alpha2,eps),@(r,y)baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi),solinit,options);

Index exceeds the number of array elements. Index must not exceed 1.

Error in solution>baseBC (line 41)
res=[ya(4); -(yRc(4))-(((m-1)/(2))*((We)^2)*((yRc(4))^3))+(yRc(6))+(((m-1)/(2))*((We)^2)*((yRc(6))^3));

Error in solution>@(r,y)baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi) (line 25)
sol = bvp4c(@(r,y)base(r,y,P0,mu,hm,Rc,We,m,Da,alpha2,eps),@(r,y)baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi),solinit,options);

Error in bvparguments (line 97)
testBC = bc(ya,yb,bcExtras{:});

Error in bvp4c (line 119)
bvparguments(solver_name,ode,bc,solinit,options,varargin);

y=deval(sol,r);
plot(r,y(1,:))

function dydx = base(r,y,P0,mu,hm,Rc,We,m,Da,alpha2,eps) % Details ODE to be solved
        %dydx = zeros(9,1);
        dydx(1)=y(4); dydx(2)=y(6); dydx(3)=y(8);
        dydx(4)=y(5); dydx(6)=y(7); dydx(8)=y(9);
        dydx(6)=(-P0)-((1/(r+eps))*(y(6)));
        dydx(8)=((-P0)/(alpha2))+((-(1)/(r+eps))*(y(8)))+((y(3))/((alpha2)*(Da)));
        dydx(4)=((1)/((1)-((mu)*(hm)*(((Rc)^2)-((r)^2)))))*(-P0+(((m-1)/(2))*((We)^2)*((3)*(y(5))*((y(4))^2)))+((mu)*(hm)*(((Rc)^2)-((r)^2))*((m-1)/(2))*((We)^2)*((3)*(y(5))*((y(4))^2)))-((2)*(mu)*(hm)*(r+eps)*(y(4)))-((2)*(mu)*(hm)*(r+eps)*((m-1)/(2))*((We)^2)*((y(4))^3))-(((1)/(r+eps))*(y(4)))-(((1)/(r+eps))*((m-1)/(2))*((We)^2)*((y(4))^3))-(((1)/(r+eps))*(mu)*(hm)*(((Rc)^2)-((r)^2))*(y(4)))-(((1)/(r+eps))*(mu)*(hm)*(((Rc)^2)-((r)^2))*((m-1)/(2))*((We)^2)*((y(4))^3)));
          % This equation is invalid at r = 0, so taking r+eps in the
          % denominator
end
    
function res = baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi) % Details boundary conditions
        res=[ya(4); -(yRc(4))-(((m-1)/(2))*((We)^2)*((yRc(4))^3))+(yRc(6))+(((m-1)/(2))*((We)^2)*((yRc(6))^3));
            yRc(1)-yRc(2); yRp(3)-yRp(2); yRp(3)-((((Da)^(1/2))/((beta)*(phi)))*((yRp(8))-(yRp(6))));
            yRb(9)-(((alpha)/((Da)^(1/2)))*((yRb(3))-((P0)*(Da)))) 
            ];
end

Answer 1

[Torsten]

@(r,y)baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi)

You don't want to use ya and yb (in your notation: r and y), namely the values of your solution variables at a and b, to define the boundary conditions ? That's impossible.

[Kalyani]

Hello Torsten.

I'm relatively new to Matlab, so forgive me if I'm not catching up with you.

In my problem, I have boundary conditions for ya, yRc, yRp, yRb, where yRb corresponds to yb (the end boundary). So I don't know what's wrong. I tried to include r and y, but I'm getting the same error again. Like I said I'm new to Matlab, so I don't know what the mistake is? Could you please elaborate?

[Torsten]

Please supply the problem you are trying to solve together with initial and boundary conditions in a mathematical formulation. Your code has so many errors that it makes no sense to correct only some settings.

[Kalyani]

Hello Torsten. I'm attaching the equations and the boundary conditions that I'm trying to solve. I've converted them to first order equations. I've provided the working in the attachment.

I'm also referring to other bvp4c codes that are available in this forum. However, I'm unable to understand where I'm going wrong. Please help.

[Kalyani]

Hello Torsten. I've tried to change r and y to ya and yRb as suggested, but I'm getting an error 'Index exceeds matrix dimensions'. I've been going through many of the bvp4c codes available in this forum and I'm not able to identify where I'm going wrong. Is the error caused by 4 boundary points instead of 2? I've seen that most of the bvp4c codes have 2 boundary points ya and yb. In my case though, I've 4 boundary points (considering 3 layers) - ya, yRc, yRp, yRb. Could this be the reason for the error? Could you give me some pointers on how to rectify the code?

[Torsten]

Is the error caused by 4 boundary points instead of 2? I've seen that most of the bvp4c codes have 2 boundary points ya and yb. In my case though, I've 4 boundary points (considering 3 layers) - ya, yRc, yRp, yRb. Could this be the reason for the error?

In this case, you have a multipoint boundary value problem. Take a look here on how to set up such a problem:

https://uk.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

[Kalyani]

Thank you,Torsten. I'll go through the link for the multipoint boundary value problem. Thank you so much!!

[Kalyani]

Hello Torsten.

I followed the link that you provided for multipoint boundary value problem. I've modified my code accordingly. The equations and the boundary conditions have been solved again and I've attached the pdf.

Here's my code:

%xc = 1;
xmesh = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
yinit = [1; 1; 1; 1; 1; 1];
sol = bvpinit(xmesh,yinit);

R=1; eps=1.0000e-04;
ya=0; yRc=0.86; yRp=0.91; yRb=1;
M=1.5;
alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
Rp=0.91;
Rb=0.95;
Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;
ud=P0*Da;

p = [R eps ya yRc yRp yRb M alpha0 P0 mu hm Rc Rp Rb Rd We m Da beta alpha alpha2 phi ud];
sol = bvp5c(@(x,y,r) f(x,y,r,p), @bc, sol);

Not enough input arguments.

Error in solution>@(x,y,r)f(x,y,r,p) (line 27)
sol = bvp5c(@(x,y,r) f(x,y,r,p), @bc, sol);

Error in bvparguments (line 96)
testODE = ode(x1,y1,odeExtras{:});

Error in bvp5c (line 126)
bvparguments(solver_name,ode,bc,solinit,options);

plot(sol.x,sol.y(1,:))

function dydx = f(x,y,region,p) % equations being solved
R = p(1);
ud = p(23);

dydx = zeros(6,1);

switch region
    case 1    % x in [0 Rc]
        dydx(1) = y(2);
        dydx(2) = (1/(1+(((m-1)/2)*((We)^2)*(3)*((y(2))^2))))*(1/(1+((mu)*(hm)*(((Rc)^2)-((x+eps)^2)))))*((-P0)+((2)*(x+eps)*(mu)*(hm)*(y(2)))+((x+eps)*(mu)*(hm)*(m-1)*((We)^2)*((y(2))^3))-((1/(x+eps))*(y(2)))-((1/(x+eps))*((m-1)/2)*((We)^2)*((y(2))^3))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*(y(2)))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*((m-1)/2)*((We)^2)*((y(2))^3)));
    case 2    % x in [Rc Rp]
        dydx(3) = y(4);
        dydx(4) = ((-P0)-((1/(x+eps))*(y(4))));
    case 3   % x in [Rp Rb]
        dydx(5) = y(6);
        dydx(6) = ((-(P0)/(alpha2))-((1/(x+eps))*(y(6)))+((y(5))/((alpha2)*(Da))));
end
end

function res = bc(YL,YR)
res = [YL(2,1)               % duc/dr=0 at r=0
       YL(2,2) + (((m-1)/2)*((We)^2)*((YL(2,2))^3)) - YL(4,2) - (((m-1)/2)*((We)^2)*((YL(4,2))^3)); %tauc=taup at r=Rc
       YL(1,2) - YL(3,2)     % uc=up at r=Rc
       YL(5,3) - YL(3,3)     % ub=up at r=Rp
       ((1/phi)*(YL(6,3))) - YL(4,3) - ((beta/(Da)^(1/2))*(YL(5,3))) %1/phi... at r=Rp
       YR(1,1) - YL(1,2)     % Continuity of uc at r=Rc
       YR(3,2) - YL(3,3)     % Continuity of up at r=Rp
       YR(5,3) - YL(5,end)   % Continuity of ub at r=Rb
       YR(6,end) - ((alpha/(Da)^(1/2))*((YR(5,end))-((P0)*(Da))))];       % dub/dr=... at r=Rb
end

I get the error 'Not enough input arguments'. Where am I going wrong?

[Torsten]

All the model parameters you define in the script part of your code are not visible in the functions for bvp5c if you don't define the functions as nested functions:

https://uk.mathworks.com/help/matlab/matlab_prog/nested-functions.html

And if you study the MATLAB example for a multipoint boundary value problem carefully, you will find that no such thing as y(3), y(4), y(5) and y(6) will exist in your case. y(1) and y(2) are just continued into the two remaining layers (possibly with different equations for them to be solved).

[Kalyani]

Hello Torsten. Thanks for your reply. I went through the link that you sent and I've modified the code, but the error persists. I get the same error 'Not enough input arguments'.

Further, you said that y(3), y(4), y(5) and y(6) will not exist. But I've second order derivatives and I'm converting them to first order. In the multipoint boundary value problem example, they have two first order odes and so only y(1) and y(2) are considered. However, for higher orders, we need to reduce to first order to use bvp4c/bvp5c. Hence, I'll get y(3), y(4), y(5) and y(6). This is based on my understanding. Please correct me if I'm wrong.

[Walter Roberson]

In my problem, I've three layers and so four boundaries - ya, yRc, yRp and yRb

Do the calculations in each layer affect the other layers (other than right at the boundary) ?

If the current position is "inside" the first layer, then are the effects of the other two layers "turned off" ? Or is it the case that the other two layers are not "turned off" but rather that they have some influence when you are fairly close to the boundary but the effects die off quickly"? Or is it more like coupled springs where each is actively affecting the others but the in-layer behavior physics is different for each layer? (Though that particular example of springs would tend to imply that the layers might change size as you proceed.)

[Kalyani]

Hello Walter.

My problem consists of equations of fluid flow in each layer and the effect of flow in one layer is felt by the nearby adjacent layer and as you said, the influence is close to the boundary.

[Walter Roberson]

If there is a point "close" to the boundary between the first and second layer, in which you cannot understand the behaviour in the first layer without knowing the current state of the second layer, then you need state variables for both layers -- possibly one state variable for the current value of each function, and one state variable for each derivative of the function except for the highest derivative. For example in the equation needs state variables for and and . Each layer might have different numbers of derivatives... keep them all .

[Kalyani]

Yes, you are right. But, I've second order terms in all the three layers. It is only on the boundaries that I need to use the state variables. I've done all these but I'm still getting error in my code.

[Walter Roberson]

Carry all of them anyhow.

Inside your ode function (and related functions), take the state vector and immediately unpack it into separate named variables, each of which are scalars

u = state(1);
du = state(2);
y1 = state(3);
dy1 = state(4);
y2 = state(5);
etc
  

Then code in terms of those variables.

n_du = du;
n_d2u = something
n_dy1 = dy1;
n_d2y1 = something;
etc

and put the output state together,

dstate = [n_du; n_d2u; n_dy1; n_d2y1 etc];

No indexing needed -- just be consistent with how you name your variables so that the meaning of the code is obvious when you read it.

[Kalyani]

Ok. But, I need the state variables only for the boundary conditions and donot need them otherwise. What about the equations? Should the equations (that are converted to first order) be kept as first order equations?

[Kalyani]

Here is my code for without using state variables.

xmesh = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
yinit = [1; 1; 1; 1; 1; 1];
sol = bvpinit(xmesh,yinit);

R=1; eps=1.0000e-04;
ya=0; yRc=0.86; yRp=0.91; yRb=1;
M=1.5;
alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
Rp=0.91;
Rb=0.95;
Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;
%ud=P0*Da;

p = [R eps ya yRc yRp yRb M alpha0 P0 mu hm Rc Rp Rb Rd We m Da beta alpha alpha2 phi];
sol = bvp5c(@(x,y,r) f(x,y,r,R,eps,ya,yRc,yRp,yRb,M,alpha0,P0,mu,hm,Rc,Rp,Rb,Rd,We,m,Da,beta,alpha,alpha2,phi), @(YL,YR)bc(YL,YR,R,eps,ya,yRc,yRp,yRb,M,alpha0,P0,mu,hm,Rc,Rp,Rb,Rd,We,m,Da,beta,alpha,alpha2,phi),sol);

plot(sol.x,sol.y(1,:))

function dydx = f(x,y,region,R,eps,ya,yRc,yRp,yRb,M,alpha0,P0,mu,hm,Rc,Rp,Rb,Rd,We,m,Da,beta,alpha,alpha2,phi) % equations being solved
R = p(1);
phi = p(22);

dydx = zeros(6,1);

switch region
    case 1    % x in [0 Rc]
        dydx(1) = y(2);
        dydx(2) = (1/(1+(((m-1)/2)*((We)^2)*(3)*((y(2))^2))))*(1/(1+((mu)*(hm)*(((Rc)^2)-((x+eps)^2)))))*((-P0)+((2)*(x+eps)*(mu)*(hm)*(y(2)))+((x+eps)*(mu)*(hm)*(m-1)*((We)^2)*((y(2))^3))-((1/(x+eps))*(y(2)))-((1/(x+eps))*((m-1)/2)*((We)^2)*((y(2))^3))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*(y(2)))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*((m-1)/2)*((We)^2)*((y(2))^3)));
    case 2    % x in [Rc Rp]
        dydx(3) = y(4);
        dydx(4) = ((-P0)-((1/(x+eps))*(y(4))));
    case 3   % x in [Rp Rb]
        dydx(5) = y(6);
        dydx(6) = ((-(P0)/(alpha2))-((1/(x+eps))*(y(6)))+((y(5))/((alpha2)*(Da))));
end
end

function res = bc(YL,YR,R,eps,ya,yRc,yRp,yRb,M,alpha0,P0,mu,hm,Rc,Rp,Rb,Rd,We,m,Da,beta,alpha,alpha2,phi)
res = [YL(2,1)               % duc/dr=0 at r=0
       YL(2,2) + (((m-1)/2)*((We)^2)*((YL(2,2))^3)) - YL(4,2) - (((m-1)/2)*((We)^2)*((YL(4,2))^3)); %tauc=taup at r=Rc
       YL(1,2) - YL(3,2)     % uc=up at r=Rc
       YL(5,3) - YL(3,3)     % ub=up at r=Rp
       ((1/phi)*(YL(6,3))) - YL(4,3) - ((beta/(Da)^(1/2))*(YL(5,3))) %1/phi... at r=Rp
       YR(1,1) - YL(1,2)     % Continuity of uc at r=Rc
       YR(3,2) - YL(3,3)     % Continuity of up at r=Rp
       YR(5,3) - YL(5,end)   % Continuity of ub at r=Rb
       YR(6,end) - ((alpha/(Da)^(1/2))*((YR(5,end))-((P0)*(Da))))];       % dub/dr=... at r=Rb
end

Now, where all the state variables need to be brought in? Because only on the boundaries, the value of the first layer and the second layer is needed. Similarly for the second and third layer. The equations in each layer, however, are distinct (they donot depend on each other). But, since the equations need to be solved with the boundary conditions and since the boundary conditions cannot be separated for each layer, I need to solve all of them simultaneously.

I'm following the example given in the multipoint boundary value problem as suggested by Torsten: https://uk.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

I've converted my second-order equations in each layer into first order equations, hence I get two first order equations in each layer. Thus, 6 first order odes need to be solved. There are 6 boundary conditions specified in the paper (converted accordingly) and three more boundary conditions are added for the continuity as given in the example of multipoint boundary value problem. There are no state variables used in the example. If I need to use state variables, how to incorporate them in the code?

[Kalyani]

I have now converted the three regions into one region ranging from 0 to 1. All the equations and boundary conditions are also converted accordingly. Here's my modified code (no different layers, since all are converted to be in a single layer ranging from 0 to 1).

function base_paper_trial_3
R=1;
ya=0; %yRc=0.86; yRp=0.91; yRb=1; 
yb=1;
%M=1.5;
%alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
Rp=0.91;
Rb=0.95;
%Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;
%ud=P0*Da;

eps=1.0000e-04;
options = [];   % place holder

sol = bvpinit(linspace(eps,1,5),[1 1 1 1 1 1]);
sol = bvp4c(@baseode,@basebc,sol,options,R,eps,ya,yb,P0,mu,hm,Rc,Rp,Rb,We,m,Da,beta,alpha,alpha2,phi);

%x = [sol.x sol.x*(lambda-1)+1];
y = [sol.y(1:2,:) sol.y(3:4,:)];

clf reset
plot(y(1,:))
%legend('v(x)','C(x)')
%title('A three-point BVP.')
%xlabel(['\lambda = ',num2str(lambda),', \kappa = ',num2str(kappa),'.'])
%ylabel('v and C')
shg

% --------------------------------------------------------------------------

function dydx = baseode(x,y,eps,P0,mu,hm,Rc,Rp,Rb,We,m,Da,alpha2)
dydx = [ Rc*(y(2))
         Rc*(1/(1+(((m-1)/2)*((We)^2)*(3)*((y(2))^2))))*(1/(1+((mu)*(hm)*(((Rc)^2)-((x+eps)^2)))))*((-P0)+((2)*(x+eps)*(mu)*(hm)*(y(2)))+((x+eps)*(mu)*(hm)*(m-1)*((We)^2)*((y(2))^3))-((1/(x+eps))*(y(2)))-((1/(x+eps))*((m-1)/2)*((We)^2)*((y(2))^3))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*(y(2)))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*((m-1)/2)*((We)^2)*((y(2))^3)));
         (Rp-Rc)*y(4)
         (Rp-Rc)*((-P0)-((1/(x+eps))*(y(4))))
         (Rb-Rp)*y(6)
         (Rb-Rp)*((-(P0)/(alpha2))-((1/(x+eps))*(y(6)))+((y(5))/((alpha2)*(Da)))) ];

% --------------------------------------------------------------------------

function res = basebc(ya,yb,P0,We,m,Da,beta,alpha,phi)
res = [ ya(2)
        (ya(2))+(((m-1)/(2))*((We)^2)*((ya(2))^3))-(ya(4))-(((m-1)/(2))*((We)^2)*((ya(4))^3))
        ya(1)-ya(3) 
        ya(3)-ya(5) 
        ((1/phi)*(ya(6)))-ya(4)-((beta/(Da)^(1/2))*(ya(5)))
        yb(6)-(((alpha)/((Da)^(1/2)))*((yb(5))-((P0)*(Da))))];

Now I get the error 'Too many input arguments'. I have removed the unwanted parameters. How to rectify this error?

[Walter Roberson]

sol = bvp4c(@baseode,@basebc,sol,options,R,eps,ya,yb,P0,mu,hm,Rc,Rp,Rb,We,m,Da,beta,alpha,alpha2,phi);

Do not use that syntax to pass "extra parameters". You do not get as much control as you would like when you use that form. I do not know if that form is even still documented -- it stopped being documented for ode45 and fmincon about 20 years ago -- and some of those functions have since added conflicting parameters so passing extra parameters that way can just fail.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Torsten]

@Kalyani

As far as I can see, you have one second-order differential equation that goes through all layers and might have a different form in each single layer. Thus you have y(1) and y(2) - the remaining four functions y(3) - y(6) do not exist because they are just y(1) and y(2) in the remaining layers.

A simple example is the heat conduction equation with different conductivities in each layer. In this case, you have two functions (T and dT/dx) in each layer, and the transmission conditions are T_left = T_right and lambda_left*dT/dx_left = lambda_right*dT/dx_right. Here, you have two functions to solve for, not 2*number_of_layers functions.

[Kalyani]

Thank you both!

You are right @Walter Roberson. I changed my code now based on the link for parametrizing functions that you sent. Thank you!

I now understand what you're saying @Torsten. I'll change my code with only y(1) and y(2) and get back to you. Thank you!!

[Kalyani]

Hello Torsten.

I've now modified my code to consist of only y(1) and y(2) as suggested. However, I'm getting the same error as before - 'Not enough input arguments'. Here's my modified code:

xmesh = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
yinit = [1; 1; 1; 1; 1; 1];
sol = bvpinit(xmesh,yinit);

R=1; eps=1.0000e-04;
ya=0; yRc=0.86; yRp=0.91; yRb=1;
M=1.5;
alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
Rp=0.91;
Rb=0.95;
Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;

p = [R eps ya yRc yRp yRb M alpha0 P0 mu hm Rc Rp Rb Rd We m Da beta alpha alpha2 phi];
sol = bvp5c(@(x,y,r) f(x,y,r,eps,P0,mu,hm,Rc,We,m,Da,alpha2), @(YL,YR)bc(YL,YR,P0,Da,beta,alpha,phi),sol);

plot(sol.x,sol.y(1,:))

function dydx = f(x,y,region,eps,P0,mu,hm,Rc,We,m,Da,alpha2) % equations being solved

dydx = zeros(6,1);
dydx(1) = y(2);
switch region
    case 1    % x in [0 Rc]
        dydx(2) = (1/(1+(((m-1)/2)*((We)^2)*(3)*((y(2))^2))))*(1/(1+((mu)*(hm)*(((Rc)^2)-((x+eps)^2)))))*((-P0)+((2)*(x+eps)*(mu)*(hm)*(y(2)))+((x+eps)*(mu)*(hm)*(m-1)*((We)^2)*((y(2))^3))-((1/(x+eps))*(y(2)))-((1/(x+eps))*((m-1)/2)*((We)^2)*((y(2))^3))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*(y(2)))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*((m-1)/2)*((We)^2)*((y(2))^3)));
    case 2    % x in [Rc Rp]
        dydx(2) = ((-P0)-((1/(x+eps))*(y(2))));
    case 3   % x in [Rp Rb]
        dydx(2) = ((-(P0)/(alpha2))-((1/(x+eps))*(y(2)))+((y(1))/((alpha2)*(Da))));
end
end

function res = bc(YL,YR,P0,Da,beta,alpha,phi)
res = [YL(2,1)               % duc/dr=0 at r=0
       YL(2,2) - YR(2,1);    % duc/dr=dup/dr at r=Rc 
       YL(1,2) - YR(1,1)     % uc=up at r=Rc
       YL(1,3) - YL(1,2)     % ub=up at r=Rp
       ((1/phi)*(YR(2,end))) - YL(2,3) - ((beta/(Da)^(1/2))*(YR(1,end))) %1/phi... at r=Rp
       YR(2,end) - ((alpha/(Da)^(1/2))*((YR(1,end))-((P0)*(Da))))];       % dub/dr=... at r=Rb
end

How to rectify the error?

[Torsten]

And note that the coordinates where one layer ends and the next layer starts have to be repeated in the grid you supply.

If you look again at the example I linked to you will see that xc appears twice in the supplied mesh:

xc = 1;
xmesh = [0 0.25 0.5 0.75 xc xc 1.25 1.5 1.75 2];

[Kalyani]

Ah..yes. I overlooked that part. Thanks for pointing it out.

But, I'm now getting another error:

Error using bvparguments (line 111)

Error in calling BVP5C(ODEFUN,BCFUN,SOLINIT):

The boundary condition function BCFUN should return a column vector of length 18.

Here's my code:

xc=0.8; xs=0.9;
xmesh = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 xc xc xs xs 1];
yinit = [1; 1; 1; 1; 1; 1];
sol = bvpinit(xmesh,yinit);

R=1; eps=1.0000e-04;
%ya=0; yRc=0.86; yRp=0.91; yRb=1;
%M=1.5;
%alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
%Rp=0.91;
%Rb=0.95;
%Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;

%p = [R eps ya yRc yRp yRb M alpha0 P0 mu hm Rc Rp Rb Rd We m Da beta alpha alpha2 phi];
sol = bvp5c(@(x,y,r) f(x,y,r,eps,P0,mu,hm,Rc,We,m,Da,alpha2), @(YL,YR)bc(YL,YR,P0,Da,beta,alpha,phi),sol);

plot(sol.x,sol.y(1,:))

function dydx = f(x,y,region,eps,P0,mu,hm,Rc,We,m,Da,alpha2) % equations being solved

dydx = zeros(6,1);
dydx(1) = y(2);
switch region
    case 1    % x in [0 Rc]
        dydx(2) = (1/(1+(((m-1)/2)*((We)^2)*(3)*((y(2))^2))))*(1/(1+((mu)*(hm)*(((Rc)^2)-((x+eps)^2)))))*((-P0)+((2)*(x+eps)*(mu)*(hm)*(y(2)))+((x+eps)*(mu)*(hm)*(m-1)*((We)^2)*((y(2))^3))-((1/(x+eps))*(y(2)))-((1/(x+eps))*((m-1)/2)*((We)^2)*((y(2))^3))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*(y(2)))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*((m-1)/2)*((We)^2)*((y(2))^3)));
    case 2    % x in [Rc Rp]
        dydx(2) = ((-P0)-((1/(x+eps))*(y(2))));
    case 3   % x in [Rp Rb]
        dydx(2) = ((-(P0)/(alpha2))-((1/(x+eps))*(y(2)))+((y(1))/((alpha2)*(Da))));
end
end

function res = bc(YL,YR,P0,Da,beta,alpha,phi)
res = [YL(2,1)               % duc/dr=0 at r=0
       YL(2,2) - YR(2,1);    % duc/dr=dup/dr at r=Rc 
       YL(1,2) - YR(1,1)     % uc=up at r=Rc
       YL(1,3) - YL(1,2)     % ub=up at r=Rp
       ((1/phi)*(YR(2,end))) - YL(2,3) - ((beta/(Da)^(1/2))*(YR(1,end))) %1/phi... at r=Rp
       YR(2,end) - ((alpha/(Da)^(1/2))*((YR(1,end))-((P0)*(Da))))];       % dub/dr=... at r=Rb
end

[Kalyani]

Now I know where the error is. Once I converted my problem in terms of only y(1) and y(2), only two initial guesses are needed and not six as before. I now get the graph. But, I don't know how to obtain a smooth curve. And, how do I obtain curves for each region separately?

Here's the corrected code:

xc=0.8; xs=0.9;
xmesh = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 xc xc xs xs 1];
yinit = [1; 1];
sol = bvpinit(xmesh,yinit);

R=1; eps=1.0000e-04;
%ya=0; yRc=0.86; yRp=0.91; yRb=1;
%M=1.5;
%alpha0=1;
P0=10;
mu=0.2;
hm=0.2;
Rc=0.86;
%Rp=0.91;
%Rb=0.95;
%Rd=1.0;
We=0.5;
m=2;
Da=0.01;
beta=0.1;
alpha=0.01;
alpha2=1.1;
phi=0.5;

%p = [R eps ya yRc yRp yRb M alpha0 P0 mu hm Rc Rp Rb Rd We m Da beta alpha alpha2 phi];
sol = bvp5c(@(x,y,r) f(x,y,r,eps,P0,mu,hm,Rc,We,m,Da,alpha2), @(YL,YR)bc(YL,YR,P0,Da,beta,alpha,phi),sol);

plot(sol.x,sol.y(1,:))

function dydx = f(x,y,region,eps,P0,mu,hm,Rc,We,m,Da,alpha2) % equations being solved

dydx = zeros(2,1);
dydx(1) = y(2);
switch region
    case 1    % x in [0 Rc]
        dydx(2) = (1/(1+(((m-1)/2)*((We)^2)*(3)*((y(2))^2))))*(1/(1+((mu)*(hm)*(((Rc)^2)-((x+eps)^2)))))*((-P0)+((2)*(x+eps)*(mu)*(hm)*(y(2)))+((x+eps)*(mu)*(hm)*(m-1)*((We)^2)*((y(2))^3))-((1/(x+eps))*(y(2)))-((1/(x+eps))*((m-1)/2)*((We)^2)*((y(2))^3))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*(y(2)))-((1/(x+eps))*(mu)*(hm)*(((Rc)^2)-((x+eps)^2))*((m-1)/2)*((We)^2)*((y(2))^3)));
    case 2    % x in [Rc Rp]
        dydx(2) = ((-P0)-((1/(x+eps))*(y(2))));
    case 3   % x in [Rp Rb]
        dydx(2) = ((-(P0)/(alpha2))-((1/(x+eps))*(y(2)))+((y(1))/((alpha2)*(Da))));
end
end

function res = bc(YL,YR,P0,Da,beta,alpha,phi)
res = [YL(2,1)               % duc/dr=0 at r=0
       YL(2,2) - YR(2,1);    % duc/dr=dup/dr at r=Rc 
       YL(1,2) - YR(1,1)     % uc=up at r=Rc
       YL(1,3) - YL(1,2)     % ub=up at r=Rp
       ((1/phi)*(YR(2,end))) - YL(2,3) - ((beta/(Da)^(1/2))*(YR(1,end))) %1/phi... at r=Rp
       YR(2,end) - ((alpha/(Da)^(1/2))*((YR(1,end))-((P0)*(Da))))];       % dub/dr=... at r=Rb
end

[Torsten]

I didn't check your conditions in detail, but I think at least

YL(1,3) - YR(1,1)     % ub=up at r=Rp

must read

YL(1,3) - YR(1,2)     % ub=up at r=Rp

And the fifth condition connects values at r = Rp with values at r = Rb. Are you sure you want this ?

YL(2,3) is at r = Rp, YR(2,end) and YR(1,end) are at r = Rb.

((1/phi)*(YR(2,end))) - YL(2,3) - ((beta/(Da)^(1/2))*(YR(1,end))) %1/phi... at r=Rp

[Kalyani]

Thank you for pointing out the error in the boundary condition and the fifth boundary condition is needed. I now get the graph like this. How do I get a smooth curve?

[Torsten]

Seems your conditions lead to a jump in the derivatives at x = 0.9. But that's what you programmed - I have no answer for your question. Check your transmission and boundary conditions if you think you must get a smooth curve throughout the complete interval. Did you correct the two conditions in your previous code ?

[Kalyani]

Did you correct the two conditions in your previous code ?

Yes, I corrected the boundary conditions and by varying some parameters, I'm able to get a smooth curve.

Thank you!!

Answer 2

[Walter Roberson]

ya is a scalar. You construct an anonymous function that ignores its parameters and passes ya to a function. Inside the called function you index ya... but ya is a scalar.

[Kalyani]

Hello Walter.

Thank you for your reply. However, could you elaborate on why ya being a scalar matters? I mean, the other boundary points yRc, yRp, yRb, are also scalars. But you are pointing out only ya as scalar which causes the error. I'm unable to understand. Could you please explain?

[Walter Roberson]

sol = bvp4c(@(r,y)base(r,y,P0,mu,hm,Rc,We,m,Da,alpha2,eps),@(r,y)baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi),solinit,options);

The second anonymous function constructed there ignores its inputs and passes the already-constructed values ya, yRc and so on to the function baseBC . At that point in the case, everything is a scalar except for r and solinit and neither of those are being passed to baseBC so every parameter to baseBC is scalar.

function res = baseBC(ya,yRc,yRp,yRb,P0,We,m,Da,beta,alpha,phi) % Details boundary conditions

The function receives the parameters under useful corresponding names -- it is not accidentally receiving a variable into a different name than you might expect. As discussed above every one of those parameters was scalar at the time the anonymous function was constructed, so each parameter received by the function will be a scalar.

        res=[ya(4); -(yRc(4))-(((m-1)/(2))*((We)^2)*((yRc(4))^3))+(yRc(6))+(((m-1)/(2))*((We)^2)*((yRc(6))^3));
            yRc(1)-yRc(2); yRp(3)-yRp(2); yRp(3)-((((Da)^(1/2))/((beta)*(phi)))*((yRp(8))-(yRp(6))));
            yRb(9)-(((alpha)/((Da)^(1/2)))*((yRb(3))-((P0)*(Da)))) 
            ];

but that code expects:

• ya to have at least 4 elements
• yRc to have at least 4 elements
• yRp to have at least 8 elements
• yRb to have at least 9 elements
• and the other parameters to be scalar

Where is ya(4) to come from when ya is a scalar that is passed unchange to the function?

[Kalyani]

Hello Walter.

Thank you for your explanation. I now understand where the problem is. But I don't know how to rectify it. Could you give me some pointers?

[Walter Roberson]

When I look through the pdf, it is not obvious to me that any indexing is needed.

I do see items such as but those appear to be functions such as so you would not be indexing them.

[Kalyani]

If I donot use indexing, then how would I define the boundary conditions? In my problem, I've three layers and so four boundaries - ya, yRc, yRp and yRb. I don't know how to define these boundaries without indexing.Can you give me some pointers?