이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Error in bvp4c code
조회 수: 22 (최근 30일)
이전 댓글 표시
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
채택된 답변
Torsten
2023년 10월 18일
@(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
Kalyani
2023년 10월 18일
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
2023년 10월 18일
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
2023년 10월 19일
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
2023년 10월 19일
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
2023년 10월 19일
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
2023년 10월 21일
Thank you,Torsten. I'll go through the link for the multipoint boundary value problem. Thank you so much!!
Kalyani
2023년 10월 22일
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
2023년 10월 22일
편집: Torsten
2023년 10월 22일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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)
Kalyani
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
편집: Torsten
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
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
2023년 10월 24일
편집: Torsten
2023년 10월 24일
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
2023년 10월 24일
편집: Torsten
2023년 10월 24일
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
2023년 10월 24일
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!!
추가 답변 (1개)
Walter Roberson
2023년 10월 18일
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
Kalyani
2023년 10월 18일
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
2023년 10월 18일
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
2023년 10월 19일
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
2023년 10월 19일
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
2023년 10월 19일
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?
참고 항목
카테고리
Help Center 및 File Exchange에서 Equation Solving에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- 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)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)