Setting up the initial conditions for BVP4C?

2 ビュー (過去 30 日間)
Richard
Richard 2014 年 5 月 25 日
編集済み: Richard 2014 年 6 月 15 日
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

回答 (1 件)

Bill Greene
Bill Greene 2014 年 5 月 28 日
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

カテゴリ

Help Center および File ExchangeLogical についてさらに検索

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by