Solving eigenvalue problem with solve_bvp

Hello everyone.

I am solving an eigenvalue problem with two unknown parameters. The colde is below:

lm0 = 8 + 11*1j;
%  Ra0 = 9000;
% options = bvpset('Stats','on','RelTol',1e-5,'NMax',10);
solinit = bvpinit(linspace(0,1,100),@mat4init,lm0);
sol = bvp4c(@mat4ode,@mat4bc,solinit);

fprintf('The fourth eigenvalue is approximately %7.3f.\n',...
        sol.parameters)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dwdx = mat4ode(z,w,lm,Ra)

chi = 1; 
m = 1.495; 
th = -chi^(-1/m)*(-1+chi^(1/m));  
l = 2;
k = 0;
a = sqrt(l^2+k^2);
% Ra = 2e05;
Pr = 0.1;
Ta = 1e05;
ph = pi/4;

eq1 = -1j*Pr*Ta^(1/2)*(l*(1+z*th)*cos(ph)-1j*m*th*sin(ph))/(Pr*(1+z*th))*w(1)-...
          Pr*Ta^(1/2)*(1+z*th)*sin(ph)*w(2)/(Pr*(1+z*th))...
          +(a^2*Pr*(1+z*th)+lm+z*th*lm)/(Pr*(1+z*th))*w(5)...
          +m*Pr*th*w(6)/(Pr*(1+z*th));
eq3 = (a^2+(1+z*th)^m*lm)*w(7)-th*w(8)/(1+z*th)-(1+z*th)^(-1+m)*w(1);
eq2 = a^2*Ra*w(7)+(2*a^2*m^2*th^2/(3*(1+z*th)^2)-(3*(m-2)*m*th^4+a^4*(1+z*th)^4)/(1+z*th)^4+lm/Pr*(-a^2-m*th^2/(1+z*th)^2)+...
      1j*k*m*Ta^(1/2)*th*cos(ph)/(1+z*th))*w(1)+(3*(m-2)*m*th^3/(1+z*th)^3+2*a^2*m*th/(1+z*th)+m*th*lm/(Pr*(1+z*th)))*w(2)+...
      (2*a^2+(4-m)*m*th^2/(1+z*th)^2+lm/Pr)*w(3)-2*m*th*w(4)/(1+z*th)+1j*l*Ta^(1/2)*cos(ph)*w(5)+Ta^(1/2)*sin(ph)*w(5);
dwdx = [ w(2)
         w(3)
         w(4)
         w(6)
         w(8)
         eq1
         eq2
         eq3];
end

function res = mat4bc(wa,wb,lm,Ra)
m = 1.495;
th = 0;
res = [  wa(1)
         wb(1)
         wa(3)+m*th/(1+th*0)*wa(2)
         wb(3)+m*th/(1+th*1)*wb(2)
         wa(6)
         wb(6)
         wa(7)
         wb(7)
         wa(2)-1
         wa(8)-1 
        ];
end

function winit = mat4init(z)
winit = [  sin(pi*z)
          0 
         0
         0
         0
         0
         0
         0];
end

But I am getting the following error

Error in mat4ode (line 19)
eq2 =
a^2*Ra*w(7)+(2*a^2*m^2*th^2/(3*(1+z*th)^2)-(3*(m-2)*m*th^4+a^4*(1+z*th)^4)/(1+z*th)^4+lm/Pr*(-a^2-m*th^2/(1+z*th)^2)+...

Error in bvparguments (line 105)
    testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 130)
    bvparguments(solver_name,ode,bc,solinit,options,varargin);
    
Error in mat4bvp (line 6)
sol = bvp4c(@mat4ode,@mat4bc,solinit);

I don't know what is wrong since the code runs well when I input Ra and calculate for lm; the problem comes when I decide to keep Ra also as an unknown parameter. Please help. Please also let me know how good convergence can be obtained and setting of initial guesses. Thank you.

[Torsten]

Take a look at this code. It shows how to properly deal with 2 unknown parameters using bvp4c.

https://de.mathworks.com/matlabcentral/answers/399320-how-to-use-bvp4c-for-boundary-value-problem-with-2-unknown-parameters

Answer 1

[Bjorn Gustavsson]

In the mat4ode function the variable Ra is undefined. How would matlab know what value to give it? What dimension should it have? 2-b-5 or 3-by-7-by-11? For the numerical evaluation all the left-hand-side variables has to be defined.

If you need to search for some special value of Ra you will have to treat the system as a function-optimization or function-root-finding problem where you solve the bvp-problem for given values of Ra and then calculates some metric/error-value/parameter of the solution that will now become a function of Ra - this function you can now use to search for the special value that otimizes/satisfies your evaluation-function.

If you can settle for a specific value of Ra - then just do that and your problem should be solved, if you change the line

sol = bvp4c(@mat4ode,@mat4bc,solinit);

to:

sol = bvp4c(@(z,w,lm) mat4ode(z,w,lm,12),@mat4bc,solinit);

You should get a solution to the problem for Ra = 12. However, the boundary-value function gives an error:

Error using bvparguments (line 111)
Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT):
  The boundary condition function BCFUN should return a column vector of length 9.
Error in bvp4c (line 130)
    bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in Some_eigenvalueproblem (line 5)
sol = bvp4c(@(z,w,lm) mat4ode(z,w,lm,12),@mat4bc,solinit); 

So you'll have to fix that.

HTH

[GAURAV]

Are you saying that bvp4c cannot be used for finding two unknown parameters?

[Bjorn Gustavsson]

This I've never tried, but if it can it will do so with an array of parameters such that you will have to merge your lm and Ra parameters into a 2-by-1 or 1-by-2 array and provide initial guesses for both and separate them in the ode- and bv-function. That is typically the way to go from 1 to multiple free parameters in matlab's optimization and ode-function.

HTH

[GAURAV]

Thank you so much for the suggestion. It is working now. However, I don't know how to get the accuracy. I am not getting the right eigenvalues right now. Also, the eigenvalues are coming out to be only real while actually they should be complex values. Below is the modified code:

lm0 = 8 + 11*1j;
Ra0 = 9000;
% options = bvpset('Stats','on','RelTol',1e-5,'NMax',10);
solinit = bvpinit(linspace(0,1,100),@mat4init,[lm0,Ra0]);
sol = bvp4c(@mat4ode,@mat4bc,solinit);
fprintf('The fourth eigenvalue is approximately %7.3f.\n',...
        sol.parameters)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dwdx = mat4ode(z,w,p)
lm = p(1);
Ra = p(2);
chi = 1; 
m = 1.495; 
th = -chi^(-1/m)*(-1+chi^(1/m));  
l = 2;
k = 0;
a = sqrt(l^2+k^2);
% Ra = 2e05;
Pr = 0.1;
Ta = 1e05;
ph = pi/4;
eq1 = -1j*Pr*Ta^(1/2)*(l*(1+z*th)*cos(ph)-1j*m*th*sin(ph))/(Pr*(1+z*th))*w(1)-...
          Pr*Ta^(1/2)*(1+z*th)*sin(ph)*w(2)/(Pr*(1+z*th))...
          +(a^2*Pr*(1+z*th)+lm+z*th*lm)/(Pr*(1+z*th))*w(5)...
          +m*Pr*th*w(6)/(Pr*(1+z*th));
eq3 = (a^2+(1+z*th)^m*lm)*w(7)-th*w(8)/(1+z*th)-(1+z*th)^(-1+m)*w(1);
eq2 = a^2*Ra*w(7)+(2*a^2*m^2*th^2/(3*(1+z*th)^2)-(3*(m-2)*m*th^4+a^4*(1+z*th)^4)/(1+z*th)^4+lm/Pr*(-a^2-m*th^2/(1+z*th)^2)+...
      1j*k*m*Ta^(1/2)*th*cos(ph)/(1+z*th))*w(1)+(3*(m-2)*m*th^3/(1+z*th)^3+2*a^2*m*th/(1+z*th)+m*th*lm/(Pr*(1+z*th)))*w(2)+...
      (2*a^2+(4-m)*m*th^2/(1+z*th)^2+lm/Pr)*w(3)-2*m*th*w(4)/(1+z*th)+1j*l*Ta^(1/2)*cos(ph)*w(5)+Ta^(1/2)*sin(ph)*w(5);
dwdx = [ w(2)
         w(3)
         w(4)
         w(6)
         w(8)
         eq1
         eq2
         eq3];
end
function res = mat4bc(wa,wb,p)
lm = p(1);
Ra = p(2);
m = 1.495;
th = 0;
res = [  wa(1)
         wb(1)
         wa(3)+m*th/(1+th*0)*wa(2)
         wb(3)+m*th/(1+th*1)*wb(2)
         wa(6)
         wb(6)
         wa(7)
         wb(7)
         wa(2)-1
         wa(8)-1 
        ];
end
function winit = mat4init(z)
winit = [  sin(pi*z)
          0 
         0
         0
         0
         0
         0
         0];
end

[Torsten]

When you call bvpinit, why do you supply only one value for lm and not a second value for Ra ?

[GAURAV]

@Torsten It was a mistake. I have corrected it.