@mmcdara 
In my code, I defined a function kk. When I substitute specific numerical values, for example kk(3, 300), the first component of the result is clearly 1.3903, which is a real number.

However, when I try to plot the function using the command:

plot(kk(3, m)[1], m = 200 .. 600)

the resulting plot is inconsistent with the expected value from the direct substitution. In particular, the y-axis values in the plot start from around 3.1, which is quite different from the 1.3903 I obtained earlier. The plotted curve doesn’t reflect the correct function values.

The same inconsistency happens when I try to create a 3D plot.

This discrepancy between direct evaluation and plotting is very confusing. I would appreciate any insights into why this is happening and how to resolve it.

@acer 

Thank you very much for your reply. However, I just tried this method, and the exported image still has a resolution of only 96 dpi, so it doesn’t seem to be a viable solution.

@mmcdara 
 

Your response is absolutely brilliant! Previously, I considered directly using the small variation at the upper endpoint to determine whether the vibration had reached a steady state. However, as you rightly pointed out, the vibration frequency may consist of multiple components, making the variation at the upper end relatively large and unreliable for identifying whether the system has entered a steady state. Using the Fourier transform is simply an excellent idea—something I had not thought of at all. I sincerely appreciate your suggestion.

That said, I’m still having some difficulty accurately identifying the steady-state region and determining the corresponding time based on this method. I gave it a try yesterday, but I haven’t yet managed to implement a frequency sweep program that incorporates the Fourier transform. Once I complete it, I’ll definitely share the results with you!

@C_R 
 

Thank you very much for your suggestion.
In previous calculations, I had indeed considered using the form of Ω = f(t) for frequency sweep analysis. For the Duffing equation, this approach successfully yielded the desired results. However, in subsequent analyses, I began to question this method. Specifically, when performing continuous calculations with sufficiently small step sizes, I became unsure about what exactly the result from the previous step represents.
As for the "steady-state" I am referring to, it mainly describes a motion state in which periodic oscillations are formed. As mentioned earlier, this corresponds to a case where the frequency content remains unchanged after performing a Fourier transform.

The current numerical analysis, as shown in the code below, implements both forward and backward frequency sweep analyses.However, when compared with the analytical solution obtained through the method of multiple scales, the numerical solution shows better agreement in the case of backward frequency sweep. For the forward sweep, the jump phenomenon tends to appear earlier than expected, which remains a challenging issue in the current analysis—quite a headache, haha. I'm also looking into the alternative method mentioned above and trying to better understand it. Once again, thank you for your detailed explanation! saopin424.mw

@C_R 
Yes, I intend to perform a frequency sweep analysis with respect to Ω based on the two differential equations mentioned above. The goal is to obtain a dataset or plot that relates Ω to the steady-state amplitude, which can then be compared with the results derived from the previously conducted method of multiple scales. The frequency sweep will be carried out over the range of angular frequency Ω = 0 to 5, during which the excitation amplitude f can be treated as a constant.

The frequency sweep analysis I envision involves using the steady-state solution at the previous frequency as the initial condition for solving at the next frequency. This approach is expected to capture the hysteresis behavior, as the solutions obtained through forward and backward sweeps may differ. If there are any misunderstandings or inappropriate assumptions in this approach, I would greatly appreciate your guidance and clarification.

@acer You wrote this code is great, it has helped me a lot, thank you very much.

@Alfred_F 
Narrow down the range of analytical solutions and then solve them in turn? That sounds pretty good,

@Alfred_F 

Yes, you are quite right that the background is indeed to solve a similar frequency and mode problem between flexible three-link

@dharr 
This is to solve the frequency of a vibration system. Although the complex solution is also part of the frequency solution, it will not be considered. Thank you very much for your reply.

@Rouben Rostamian  
Thank you for your answer, I didn't think of using the form of graph combination to solve, thank you very much, and found that the solution of fsolve is much faster than analytic

@acer  This equation should have several non-trivial solutions, but I'm not sure about the initial range

@dharr  Ok, I'm going to try it out

@acer 

I used maple to compute the solution of a polynomial equation, which takes a long time and is slow to draw graphs, using the odeplot command at that time

@dharr I've solved it in a different way, I've removed one of the variables, I've got a new solution, I don't know how to solve it, or is there something wrong with my equation.

question11.mw

@dharr 

Why doesn't solve make sense? I was hopeful when I faced this result, but is there no nontrivial solution to this equation? Can we get rid of the trigonometric function and then solve it?