Third-order polynomial equation which complex roots

Hi,
I want to plot the three roots of c (real and imaginary) as a function of k for the following third-order polynomial equation:
I am using fsolve to code it but this requires three initial guesses which are hard to identify for the given equation.
Any suggestions?
Thank you

 Respuesta aceptada

Carola Forlini
Carola Forlini el 6 de Feb. de 2024
Thank you for all the answers.
At the end the easiest way was to calculate first the discriminant of thee polynomio for a range of k and then use the roots function to calculate the solution. In this way I have better control on the expected solutions since the discriminant will tell me if I should have all real roots or real and complex conjugates one.

Más respuestas (2)

Dyuman Joshi
Dyuman Joshi el 3 de Feb. de 2024
Movida: Matt J el 3 de Feb. de 2024
Define the polynomial as a function handle of the variable 'k' and use roots for different values of 'k'.
Also, note that you will need to plot the real and imaginary separately.
syms c k L
eqn = c^3 ...
- c^2*(2*k*exp(4*k*L) + 2*k*exp(2*k*L) + exp(4*k*L) - 6*k^2*exp(2*k*L) + 1) / (2*exp(2*k*L) * k * (exp(2*k*L)+1)) ...
- c*(-k*exp(4*k*L) + 2*k*exp(2*k*L) - k + 2*exp(4*k*L) - 2) / (2*exp(2*k*L)*k^2*(exp(2*k*L) + 1)) ...
+ (exp(4*k*L) - 2*exp(2*k*L) + 1) / (2*exp(2*k*L)*k^3*(exp(2*k*L) + 1))
eqn = 
solutions = solve(eqn, c, 'maxdegree', 3)
solutions = 
sol= subs(solutions, L, 2); %arbitrary
%vpa(limit(sol(1), k, 0, 'left'))
%vpa(limit(sol(1), k, 0, 'right'))
tiledlayout('flow');
nexttile(); fplot([real(sol(1)), imag(sol(1))], [-3 3]); title('root #1');
nexttile(); fplot([real(sol(2)), imag(sol(3))], [-5 5]); title('root #2');
nexttile(); fplot([real(sol(3)), imag(sol(3))], [-3 3]); title('root #3');

3 comentarios

syms c k L
eqn = c^3 ...
- c^2*(2*k*exp(4*k*L) + 2*k*exp(2*k*L) + exp(4*k*L) - 6*k^2*exp(2*k*L) + 1) / (2*exp(2*k*L) * k * (exp(2*k*L)+1)) ...
- c*(-k*exp(4*k*L) + 2*k*exp(2*k*L) - k + 2*exp(4*k*L) - 2) / (2*exp(2*k*L)*k^2*(exp(2*k*L) + 1)) ...
+ (exp(4*k*L) - 2*exp(2*k*L) + 1) / (2*exp(2*k*L)*k^3*(exp(2*k*L) + 1))
eqn = 
solutions = solve(eqn, c, 'maxdegree', 3)
solutions = 
sol= subs(solutions, L, 2); %arbitrary
%vpa(limit(sol(1), k, 0, 'left'))
%vpa(limit(sol(1), k, 0, 'right'))
tiledlayout('flow');
nexttile(); fplot([real(sol(1)), imag(sol(1))], [-3 -1e-2]); title('root #1');
nexttile(); fplot([real(sol(1)), imag(sol(1))], [1e-2, 3]); title('root #1');
%nexttile(); fplot([real(sol(2)), imag(sol(3))], [-5 5]); title('root #2');
%nexttile(); fplot([real(sol(3)), imag(sol(3))], [-3 3]); title('root #3');
syms c k L
eqn = c^3 ...
- c^2*(2*k*exp(4*k*L) + 2*k*exp(2*k*L) + exp(4*k*L) - 6*k^2*exp(2*k*L) + 1) / (2*exp(2*k*L) * k * (exp(2*k*L)+1)) ...
- c*(-k*exp(4*k*L) + 2*k*exp(2*k*L) - k + 2*exp(4*k*L) - 2) / (2*exp(2*k*L)*k^2*(exp(2*k*L) + 1)) ...
+ (exp(4*k*L) - 2*exp(2*k*L) + 1) / (2*exp(2*k*L)*k^3*(exp(2*k*L) + 1));
solutions = solve(eqn, c, 'maxdegree', 3);
sol= subs(solutions, L, 2); %arbitrary
%vpa(limit(sol(1), k, 0, 'left'))
%vpa(limit(sol(1), k, 0, 'right'))
tiledlayout('flow');
%nexttile(); fplot([real(sol(1)), imag(sol(1))], [-3 3]); title('root #1');
%nexttile(); fplot([real(sol(2)), imag(sol(3))], [-5 5]); title('root #2');
nexttile(); fplot([real(sol(3)), imag(sol(3))], [-1 1]); title('root #3');
syms c k L
eqn = c^3 ...
- c^2*(2*k*exp(4*k*L) + 2*k*exp(2*k*L) + exp(4*k*L) - 6*k^2*exp(2*k*L) + 1) / (2*exp(2*k*L) * k * (exp(2*k*L)+1)) ...
- c*(-k*exp(4*k*L) + 2*k*exp(2*k*L) - k + 2*exp(4*k*L) - 2) / (2*exp(2*k*L)*k^2*(exp(2*k*L) + 1)) ...
+ (exp(4*k*L) - 2*exp(2*k*L) + 1) / (2*exp(2*k*L)*k^3*(exp(2*k*L) + 1));
solutions = solve(eqn, c, 'maxdegree', 3);
sol= subs(solutions, L, 2); %arbitrary
vpa(limit(sol(3), k, 0, 'left'))
ans = 
vpa(limit(sol(3), k, 0, 'right'))
ans = 

Iniciar sesión para comentar.

Categorías

Más información sobre Symbolic Math Toolbox en Centro de ayuda y File Exchange.

Preguntada:

el 3 de Feb. de 2024

Respondida:

el 6 de Feb. de 2024

Community Treasure Hunt

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

Start Hunting!

Translated by