Error in bvp4c code

Question

Kalyani 2023-10-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2035139-error-in-bvp4c-code

评论： Kalyani 2023-10-24

采纳的回答： Torsten

equations_boundary_conditions.pdf

在 MATLAB Online 中打开

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

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

请先登录，再进行评论。

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

Answer 1

Torsten 2023-10-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2035139-error-in-bvp4c-code#answer_1335774

在 MATLAB Online 中打开

@(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 个评论
显示 26更早的评论隐藏 26更早的评论

Kalyani 2023-10-22

在 MATLAB Online 中打开

base_equations_boundary_conditions.pdf

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?

Kalyani 2023-10-24

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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?

Kalyani 2023-10-24

在 MATLAB Online 中打开

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?

Kalyani 2023-10-24

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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

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

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!!

请先登录，再进行评论。

Answer 2

Walter Roberson 2023-10-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2035139-error-in-bvp4c-code#answer_1335564

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 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2023-10-18

在 MATLAB Online 中打开

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 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.

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?

请先登录，再进行评论。

Error in bvp4c code

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

采纳的回答

28 个评论
显示 26更早的评论隐藏 26更早的评论

更多回答（1 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Error in bvp4c code

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

采纳的回答

28 个评论 显示 26更早的评论隐藏 26更早的评论

更多回答（1 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

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

28 个评论
显示 26更早的评论隐藏 26更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论