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

7 views (last 30 days)
Tanya Sharma on 3 Oct 2019
Commented: darova on 11 Oct 2019
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)
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?

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);
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
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

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.
darova on 11 Oct 2019
Cool. It works!

R2018b

### Community Treasure Hunt

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

Start Hunting!

Translated by