1
$\begingroup$

I want to solve the Duffing equation, which is a second-order differential equation. I am able to get the steady-state plot, but I do not know how to get maximum value from the last few points and the vary the frequency to get a plot of w vs y. I have tried the following code to just startup

ClearAll[m, k, F, w, k3, e, g]
m = 1;
k = 1;
F = 0.01;
w = 1;
k3 = 10^-4;
e = 10^-5;
g = 0.001;
center = 1;
range = .05;
start = center - range/2;
stop = center + range/2;
pts = 100;
step = (stop - start)/pts;
eqn = {m*y''[x] + m*g*y'[x] + k*y[x] + k3*y[x]^3 + e*y[x]^2*y'[x] - 
     F*Cos[w*x] == 0, y[0] == 0, y'[0] == 0};
c = Range[start, stop, step];
s = Table[NDSolve[eqn, y, {x, 0, 50}], {w, Length@c}];
Plot[s, {t, 0, 60}, PlotLegends -> "Expressions", ImageSize -> Large]
$\endgroup$
1
  • 2
    $\begingroup$ Do you know about ParametricNDSolve? I think it does what you need. $\endgroup$
    – Roman
    Nov 30, 2020 at 14:44

1 Answer 1

0
$\begingroup$

While the comment of Roman (using ParametricNDSolve) is definitely a valid point and most likely the best solution for this kind of application. The manual approach suggested by OP is not impossible. The code presented by OP has several typos in the line defining s, which I fixed in the version below

ClearAll[m,k,F,w,k3,e,g]
m=1;
k=1;
F=0.01;
k3=10^-4;
e=10^-5;
g=0.001;
center=1;
range=.05;
start=center-range/2;
stop=center+range/2;
pts=100;
step=(stop-start)/pts;
eqn={m*y''[x]+m*g*y'[x]+k*y[x]+k3*y[x]^3+e*y[x]^2*y'[x]-F*Cos[w*x]==0,y[0]==0,y'[0]==0};
c=Range[start,stop,step];
s=Table[NDSolveValue[eqn,y,{x,0,50}][t],{w,c}];
p1=Plot[s,{t,0,50},ImageSize->Large,PlotRange->All,FrameLabel->{"x","y[x]"},GridLines->Automatic];

wData=MapThread[{#1,NMaximize[{#2,0<=t<=50},t]}&,{c,s}];
Transpose[{wData[[All,1]],wData[[All,2,1]]}];
p2=ListLinePlot[%,Frame->True,FrameLabel->{"w","Subscript[y, max]"},GridLines->Automatic,ImageSize->Large];
SortBy[%%,-#[[2]]&]//First
Grid[{{p1,p2}}]

which results in (Warning: the computation took around 40 seconds of wall time on my laptop)

Output

where the maximum amplitude of $\sim 0.2259$ is reached for $w=0.9945$. The same could be realized using ParametricNDSolve in the code above:

...
eqn={m*y''[x]+m*g*y'[x]+k*y[x]+k3*y[x]^3+e*y[x]^2*y'[x]-F*Cos[w*x]==0,y[0]==0,y'[0]==0};
c=Range[start,stop,step];
sol=ParametricNDSolveValue[eqn,y,{x,0,50},w]
s=Table[sol[w][t],{w,c}];
...

which took again roughly 40 seconds to evaluate (I have not benchmarked this using multiple runs). Which makes sense since ParametricNDSolveValue and its output ParametricFunction is basically a fancy wrapper for NDSolve (with maybe some additional functionalities).

$\endgroup$
3
  • $\begingroup$ @Nova, please correct me if I am wrong. what is the role of t in this line s=Table[NDSolveValue[eqn,y,{x,0,50}][t],{w,c}]; $\endgroup$ Dec 1, 2020 at 7:31
  • $\begingroup$ @ShelenderKumar The line computes the solutions for the given set of w: stored in c up to t=50. One can only reliably plot/evaluated the computed solutions between t=0 and t=50. If one wants to evolve to later times one would need to change {x,0,50} in NDSolveValue to something like {x,0,tmax}. $\endgroup$
    – N0va
    Dec 1, 2020 at 13:04
  • $\begingroup$ I want to change the initial condition in such a way that every time a new value of c is taken initial condition changes to y(0)=y(end) and y'(0)=y'(end);. $\endgroup$ Dec 1, 2020 at 17:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.