duffing_lindstedt

As always, start with resetting the workspace.

reset()

The original differential equation.

diffeq:=diff(x(t),t,t)+x(t)+epsi*x(t)^3

diff(x(t), t, t) + epsi*x(t)^3 + x(t)

We do a transformation such that tau = omega*t

diffeq:=omega^2*diff(x(tau),tau,tau)+x(tau)+epsi*x(tau)^3

omega^2*diff(x(tau), tau, tau) + epsi*x(tau)^3 + x(tau)

omega is the frequency, which has an epsilon dependence.

omega:=1+k1*epsi

x(tau):=x0(tau)+epsi*x1(tau)

diffeq:=Simplify(diffeq)

x0(tau) + epsi*(x0(tau) + epsi*x1(tau))^3 + epsi*x1(tau) + (epsi*k1 + 1)^2*(epsi*diff(x1(tau), tau, tau) + diff(x0(tau), tau, tau))

diffeq_0:=coeff(diffeq,epsi,0)

diff(x0(tau), tau, tau) + x0(tau)

diffeq_1:=coeff(diffeq,epsi,1)

x1(tau) + 2*k1*diff(x0(tau), tau, tau) + x0(tau)^3 + diff(x1(tau), tau, tau)

diffeq_0:=ode({diffeq_0=0,x0(0)=A,x0'(0)=0},x0(tau))

x0_solution:=solve(diffeq_0); x0_solution:=x0_solution[1];

diffeq_1:=Simplify(subs(diffeq_1,x0(tau)=x0_solution))

diff(x1(tau), tau, tau) + cos(tau)*(cos(2*tau)/2 + 1/2)*A^3 - 2*k1*cos(tau)*A + x1(tau)

diffeq_1:=ode({diffeq_1=0,x1'(0)=0,x1(0)=0},x1(tau))

$ode({D(x1)(0) = 0, x1(0) = 0, (D@@2)(x1)(tau) + x1(tau) + cos(tau)*(cos(2*tau)/2 + 1/2)*A^3 - 2*k1*cos(tau)*A}, x1(tau))$

x1_solution:=solve(diffeq_1)

${- cos(tau)*((A*k1)/2 + (A^3*cos(2*tau))/8 + (A^3*cos(4*tau))/32 - (5*A^3)/32 - (A*k1*cos(2*tau))/2) - sin(tau)*((A^3*sin(2*tau))/4 + (A^3*sin(4*tau))/32 + (3*A^3*tau)/8 - (A*k1*sin(2*tau))/2 - A*k1*tau)}$

x1_solution:=Simplify(x1_solution[1])

-(A*sin(tau)*((A^2*sin(2*tau))/2 - 8*k1*tau + 3*A^2*tau))/8

collect(x1_solution,tau)

((A*sin(tau)*(8*k1 - 3*A^2))/8)*tau - (A^3*sin(2*tau)*sin(tau))/16

The terms in the bracket, multiplying tau, consist of the secular term. We would like it to be zero. And now, we have the freedom to make it zero.

secularterm:=coeff(x1_solution,tau,1)

A*k1*sin(tau) - (3*A^3*sin(tau))/8

We add a bunch of assumptions about the various variables. Try deleting them one by one to see what their effect of k1 is. MuPad will find more solutions assuming every variable to be generally complex-valued.

assume(k1 in R_); assume(tau in R_); assume(A in R_); assumeAlso(sin(tau) <> 0); assumeAlso(A <>0);

k1:=solve(secularterm,k1);

${(3*A^2)/8}$

k1:=k1[1]

(3*A^2)/8

Now let us look at the frequency 'omega'. By using the above k1 in omega, the following equation gives the relationship between the natural frequency 'omega' and the amplitude A.

omega

(3*epsi*A^2)/8 + 1

As we see above, for epsi>0, the frequency omega increases with amplitude. For epsi<0 (as would be the case if we approximated the simple pendulum), the frequency omega decreases and the time period increases with amplitude. This agrees with what we know about the simple pendulum.

x1_solution

-(A^3*sin(2*tau)*sin(tau))/16

x_solution:=x0_solution+epsi*x1_solution

A*cos(tau) - (A^3*epsi*sin(2*tau)*sin(tau))/16

tau:=omega*t

t*((3*epsi*A^2)/8 + 1)

x_solution

A*cos(t*((3*epsi*A^2)/8 + 1)) - (A^3*epsi*sin(t*((3*epsi*A^2)/8 + 1))*sin(2*t*((3*epsi*A^2)/8 + 1)))/16

xdot_solution:=diff(x_solution,t)

- A*sin(t*((3*epsi*A^2)/8 + 1))*((3*epsi*A^2)/8 + 1) - (A^3*epsi*sin(t*((3*epsi*A^2)/8 + 1))*cos(2*t*((3*epsi*A^2)/8 + 1))*((3*epsi*A^2)/4 + 2))/16 - (A^3*epsi*sin(2*t*((3*epsi*A^2)/8 + 1))*cos(t*((3*epsi*A^2)/8 + 1))*((3*epsi*A^2)/8 + 1))/16

epsi:=0.1;

for i from 1 to 10 step 1 do
A:=i*0.1;
temp1:=subs([x_solution,xdot_solution],A=i*0.1);
curve1[i]:=plot::Curve2d(temp1, t=0..(2*PI/omega));
end_for

plot(curve1[i]$i=1..10)

MuPAD graphics

x_solution