Setting up the initial conditions for BVP4C?

3 次查看（过去 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 Exchange 中查找有关 MATLAB Mobile Fundamentals 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by