I'm trying to use BVP4C in a loop where the answer from the previous loop is the initial conditions for the next loop. Do I still use the Solinit function to get the initial conditions in the right format or should I do something else? Also, once the code gets to BVP4C, I get the following error:

Error using bvparguments (line 109) Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT,P1,P2...): The derivative function ODEFUN should return a column vector of length 2.

I've seen some examples where the vector separator is a semicolon, but others where it appears to be a space. What exactly should I use?

The for loop that I used for testing is shown below:

%initial conditions

p(1) = pratio;

for i = 2:pnts

p(i)= ((1/pnts/10)*i) + pratio; %using x dimension

Ma(i) = 0.0;

end

%start of full calculation loop, start with 2 iterations

%while (err>1*10^(-5))

for j = 1:2

ptot=0.0;

x=linspace(0,1,pnts);

%calculating q(x) eq. 26

for i = 1:pnts

Ao = 0.0;

A_n =0.0;

eta=acos(i/pnts);

cons = 1.0/sqrt((sinh(eps1)*sinh(eps1))+sin(eta)*sin(eta));

for k =1:200

An = 0.0;

%to calculate An and Ao

%calculating the integral along fracture

for m = 1:(pnts-1)

test =cos(2*acos(m/pnts));

An = An + (1.0 - (p(m)*p(m)))*(cos(2*k*(acos(m/pnts)))/...

sqrt(1.0 - m*m/pnts/pnts));

if k==1

Ao = Ao + (1.0 - (p(m)*p(m)))*(1.0 - m*m/pnts/pnts)^(-0.5);

else

end

end

An = An*4.0/pi/(sinh(2.0*k*(eps_inf - eps1)));

A_n = A_n + An;

A_nn(k) = An;

end

Ao = (-2.0)/pi/(eps_inf - eps1)*Ao;

q(i)= lamda*cons/p(i)*(A_n - Ao);

end

plot(x,q)

% For fixed well pressure:

%param(6)=pratio; % if used for fixed well pressure condition

% How do I get q(x) in the right form?

solinit=bvpinit(linspace(0,1,pnts),[0.01; 0.9]); %0.01 intial Mach #, 0.9 initial pressure

%C3= (1.0 - gama*(y(1)^2.0)); % C3

sol=bvp4c(@frac,@fracbc,solinit);

%x=linspace(0,1,pnts);

y=deval(sol,x);

%store last values for error calculation

for i = 1:pnts

p_1(i)= p(i);

Ma_1(i) = Ma(i);

ptot=ptot+p(i);

end

%write calculated values to arrays

for i = 1:pnts

p(i)= y(2,i);

Ma(i) = y(1,i);

end

%end of the while loop, for testing: a for loop

end

Here are the two nested functions that I forgot to put in:

function res=fracbc(ya,yb) %#ok

% y(1) is Mach #

% y(2) is p, pressure

Ao = 0.0;

A_n =0.0;

%calculating each the infinite series to 200 terms

for k =1:200

An = 0.0;

%to calculate An and Ao

%calculating the integral along fracture

for m = 1:(pnts-1)

An = An + (1.0 - (p(m)*p(m)))*(cos(2*k*(acos(m/pnts)))/...

sqrt(1.0 - m*m/pnts/pnts));

if k==1

Ao = Ao + (1.0 - (p(m)*p(m)))*(1.0 - m*m/pnts/pnts)^(-0.5);

else

end

end

An = An*4.0*k*(cosh(2.0*k*(eps_inf - eps1)))...

/pi/(sinh(2.0*k*(eps_inf - eps1)));

A_n = A_n + An;

A_nn(k) = An;

end

Ao = (-2.0)/pi/(eps_inf - eps1)*Ao;

M = (1.0/alpha/cfd/p(pnts))*(2.0*A_n - Ao);

res = [ya(2)-pratio;yb(1)-M]; % This is for fixed pressure at well

end

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

function dydx=frac(x,y)

C3 = (1.0 - gama*(y(1).^2.0));

C2 =sqrt(gama);

M_P = y(1)./y(2);

B = (-y(1).*alpha - beta*y(2).*y(1).*y(1) - C2*q.*y(1).*y(2))/C3;

A = (q./C2) - M_P.*B;

dydx=[A;B];

end

Answer by Bill Greene
on 28 May 2014

From your call to bvpinit where you have this [0.01; 0.9] for the initial guess, it looks like you have 2 differential equations in your system. (In MATLAB, the ; means the entries that follow are in the next row so this is a column vector of length two.)

You didn't include the code for your frac function but apparently it isn't returning a column vector of length two, judging from the error message.

In MATLAB, there are two common ways to construct column vectors:

[1; 2; 3; 4] % create a column vector of length four directly

or

[1 2 3 4]' % create a row vector and then transpose it

Bill

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.