Solving eigenvalue problem with solve_bvp
3 views (last 30 days)
Show older comments
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.
1 Comment
Torsten
on 11 Jan 2022
Take a look at this code. It shows how to properly deal with 2 unknown parameters using bvp4c.
Answers (1)
Bjorn Gustavsson
on 11 Jan 2022
Edited: Bjorn Gustavsson
on 11 Jan 2022
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
5 Comments
Torsten
on 11 Jan 2022
When you call bvpinit, why do you supply only one value for lm and not a second value for Ra ?
See Also
Categories
Find more on Boundary Value Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!