How can I solve a system of ODEs having coefficients in vector form using bvp4c ?

3 views (last 30 days)
Hi,
I am solving a system of 5 ODEs that contains known parameters, unknown parameters and some coefficients are of the form of vector.
%-------------------------------ODE system-----------------------------------%
function eqns = odes(x,y,e)
global Pr phi Ra Da Fr A1 A2 fdesh fdeshdesh thetadesh;
eqns = [y(2)
y(3)
(phi./Da).*y(2)+(2.*phi.*Fr./A1).*fdesh.*y(2)-(fdesh.*1./A1).*y(3)-(fdeshdesh.*1./A1).*y(1)+(2.*fdesh.*1./A1).*y(2)-(e./A1).*y(2)-(phi.*Ra./(A1^2).*A2).*y(4)
y(5)
-(Pr./A2).*(fdesh.*y(5)+thetadesh.*y(1)+e.*y(4))];
end
%------------------------------------------------------------------------------------------------
--------------------% fdesh, fdeshdesh and thetadesh are of vector form. these are the known solution of the governing equations .%----------------------
%----------------------------------------------------------------------------------------------
%--------------------------boundary conditions-----------------------------%
function res = ode_bc(ya,yb,e)
res = [ya(1)
ya(2)
ya(3)
ya(4)
yb(2)
yb(4)];
end
%-------------------------------------------------------------------------------------------------
I am getting this error:
Error using bvparguments (line 108)
Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT):
The derivative function ODEFUN should return a column vector of length 5.
%-------------------------------------------------------------------------------------------------------
are the known solutions causing a problem to solve the system?
Thanks in advance.

Accepted Answer

darova
darova on 3 Oct 2019
Try this
function eqns = odes(x,y,e,x0)
global Pr phi Ra Da Fr A1 A2 fdesh fdeshdesh thetadesh; % global variables are not recommended
% x0 - vector of corresponding values for fdesh, fdeshdesh, thetadesh
% x0(end) should not be bigger than x(end)
fd = interp1(x0,fdesh,x);
fdd = interp1(x0,fdeshdesh,x);
thd = interp1(x0,thetadesh,x);
eqns = [y(2)
y(3)
(phi./Da).*y(2)+(2.*phi.*Fr./A1).*fd.*y(2)-(fd.*1./A1).*y(3)-(fdd.*1./A1).*y(1)+(2.*fd.*1./A1).*y(2)-(e./A1).*y(2)-(phi.*Ra./(A1^2).*A2).*y(4)
y(5)
-(Pr./A2).*(fd.*y(5)+thd.*y(1)+e.*y(4))];
end
  10 Comments
Tanya Sharma
Tanya Sharma on 9 Oct 2019
Indeed it is darova.
Still want to find the unkowns "e".
What if the whole problem is a linearized eigenvalue problem, where the "e" are the unknown eigenvalues. In that case would I be able to trace the eigenvalues "e" using bvp4c ?
darova
darova on 9 Oct 2019
I don't think β patameter can be found. In the paper you attached everywhere is said that it can be obtained with guess

Sign in to comment.

More Answers (1)

Tanya Sharma
Tanya Sharma on 10 Oct 2019
Yes and the guess is provided in the "solinit" structure as an extra parameter
beta = 99;
init2 = bvpinit(linspace(-1,1,50),@math4init,beta);
Which we can later be obtain by "sol.parameters", this provides us the value bvp4c has taken nearer to our provided guess.
Which in this case is:
beta=77.8263.
  3 Comments
Tanya Sharma
Tanya Sharma on 11 Oct 2019
Hi,
I did a little change and it is passing correctly
sol = bvp4c(@(x,y,e)odes(x,y,data,params,e),@(ya,yb,e)ode_bc(ya,yb,e), solinit);
fprintf('e=%g.\n',sol.parameters);
which returns:
e=1.00092.
But getting a warning:
Warning: Unable to meet the tolerance without using more than 2000 mesh points.
The last mesh of 2500 points and the solution are available in the output argument.
The maximum residual is 0.00567401, while requested accuracy is 0.001.
> In bvp4c (line 323)
In exer1 (line 76)
May be have to change the initial guess or th bcs. Otherwise its fine.

Sign in to comment.

Categories

Find more on Parallel Computing Fundamentals in Help Center and File Exchange

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!