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

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.

Answer 1

[darova]

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

[Tanya Sharma]

Thanks darova.

I am still having some issues with the code.As it is giving me intepolation error on interp1. saying not enough input arguments

1. as you suggested that x0 must be the vector of corresponding values for fdesh,sdeshdesh,thetadesh. But the 'fdesh,sdeshdesh,thetadesh' out of the ode function are being evaluated at x. that i have taken linspace(1,20,100). So does that mean it should be defined as fd=interp1(x,fdesh,x0). where x0 must be coarser than x?
2. Also you did not recommend the global variables. So how would I define the all the parameters and the known fdesh,sdeshdesh,thetadesh inside the odes function?

Thanks again.

[darova]

You don't know how x vector dense/coarse is (it is inside function). Vector x that ode45() returns can be any step (try to change it outside ode45, it doesn't influence on precision)

fd = interp1(x0,fdesh,x);

To pass variable or array you use something like:

data = [x0(:) fdesh(:) fdeshfdesh(:) thetadesh(:)];
params = [Pr phi Ra Da Fr A1 A2 e];
%%
function eqns = odes(x,y,data,params)
% global Pr phi Ra Da Fr A1 A2 fdesh fdeshdesh thetadesh;        % global variables are not recommended
x0 = data(:,1);
%% ---
thetadesh = data(:,4);
Pr = params(1);
%% ---
e = params(8);

If those parameters (Pr phi Ra Da Fr A1 A2) are arrays too, then you should use interp1()

[Tanya Sharma]

Thanks again drova,

I am following the procedure you gave me above but still getting some problems while parsing the values of "data" and "params" inside the odes( ) function:

1. error is occuring in parsing the data(:,2). As it says "Index in position 2 exceeds array bounds (must not exceed 1).". But otherwise it is working fine while I am using it in the script.
2. I want to convey that I am using the bvp4c solver not ode45.
3. I am also still confused about the mesh that I need to use, x and x0. As the f,fdesh,fdeshdesh,thetadesh in the script are evaluated using the 'x' mesh. So won't it needed to be interpolated using a mesh x0 that may be witten as "fd=interp1(x,fdesh,x0);" I am not sure if I am right or not but still a beginner.

Tanya.

[darova]

You are almost there!

There is something wrong with your boundary condition (there are too much of them)

Also i changed initial conditions (because all solutions were zeros)

See attached script

[Tanya Sharma]

Thanks a lot again darova,

The code is running well with this. I am resolved for the half of my problem with this.

The next thing is :

1. I want to evaluate the values of unknown parameter "e" in the equations, which is considered as the eigenvalue.
2. The boundary conditions that I supplied was an extra condition, because of the presence of an unknown parameter "e" in the equations. Do we need to provide an extra condition for the unkown parameter(eigenvalue) "e"?

[darova]

Do we need to provide an extra condition for the unkown parameter(eigenvalue) "e"?

You tell me! What is condition for 'e' parameter?

[Tanya Sharma]

See page 180 of the attached pdf.

I am supplying the extra condition "ya(3)".

[darova]

Beauty?

function res = ode_bc(ya,yb,e)
    res = [ya(1)
           ya(2)
           ya(3)-0.1
%            ya(4)
           yb(2)
           yb(4)];
end

[Tanya Sharma]

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]

I don't think β patameter can be found. In the paper you attached everywhere is said that it can be obtained with guess

Answer 2

[Tanya Sharma]

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]

How to pass it correctly?

[Tanya Sharma]

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.

[darova]

Cool. It works!