You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Error in bvp4c code
22 views (last 30 days)
Show older comments
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
Accepted Answer
Torsten
on 18 Oct 2023
@(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.
28 Comments
Kalyani
on 18 Oct 2023
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?
Kalyani
on 19 Oct 2023
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
on 19 Oct 2023
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
on 19 Oct 2023
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:
Kalyani
on 22 Oct 2023
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
on 22 Oct 2023
Edited: Torsten
on 22 Oct 2023
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:
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
on 24 Oct 2023
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
on 24 Oct 2023
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.)
Walter Roberson
on 24 Oct 2023
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 .
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1519726/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1519731/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1519736/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1519741/image.png)
Walter Roberson
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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.
Torsten
on 24 Oct 2023
Edited: Torsten
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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
on 24 Oct 2023
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,:))
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1520071/image.png)
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
on 24 Oct 2023
Edited: Torsten
on 24 Oct 2023
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
Torsten
on 24 Oct 2023
Edited: Torsten
on 24 Oct 2023
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 ?
More Answers (1)
Walter Roberson
on 18 Oct 2023
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.
5 Comments
Kalyani
on 18 Oct 2023
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
on 18 Oct 2023
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?
Walter Roberson
on 19 Oct 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1515059/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1515064/image.png)
Kalyani
on 19 Oct 2023
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?
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)