Plotting first and second order ODE

726 Views Asked by At

Here is the second order differential equation:

$\frac{d^2x}{dt^2} + \sin(x) = 0$

with initial conditions $x(0)=1$ and $x'(0)=0$.

The written code in matlab that numerically solves and plots the results is shown below:

function second_order
t=0:0.001:14;
initial_x=1;
initial_dxdt = 0;
[t,x]=ode45(@rhs,t,[initial_x initial_dxdt]);
plot(t,x(:,1));
xlabel('t');ylabel('x');

function dxdt=rhs(t,x)
dxdt_1=x(2);
dxdt_2=-sin(x(1));
dxdt = [dxdt_1; dxdt_2);
end
end

The plot shows the sinusoidal curve with maximum amplitude plus and minus one. This plot is correct. However, when I turn the second order differential equation into first one as follows:

$ \frac{d^2x}{dt^2}\frac{dx}{dt} + \sin(x)\frac{dx}{dt} = 0$

$\frac{1}{2}(\frac{dx}{dt})^2 = \cos(x)-\cos(1)$ and that leads to first order ODE,

$\frac{dx}{dt}=\sqrt{2(\cos(x)-\cos(1))}$

The matlab to plot first ODE is shown below:

function first_order
t=0:0.001:14;
initial_x=1;
[t,x]=ode45(@rhs, t, initial_x);
plot(t,x);
xlabel('t'); ylabel('x');

function dxdt=rhs(t,x)
dxdt = sqrt(2)*sqrt(cos(x)-cos(1));
end
end

The problem is that I am not getting the same plot as in second ODE code. Instead, I am getting a straight line at $x=1$. I know that the code is correct because I tested it with other first order differential equation. Therefore, why I am not getting the same plot even thought the first and second order differential equations are the same. First order is basically a solution of the second order. But values of $x$ should be the same. Is this approach not applicable? or am I doing something wrong?

3

There are 3 best solutions below

1
On BEST ANSWER

You can get the "curved" solution from a first order equation by further parametrizing the integrated equation $$ \frac12\dot x^2+\frac12\left(2\sin\frac x2\right)^2=\frac12\left(2\sin\frac 12\right)^2 $$ by recognizing this as a circle equation for $(\dot x,2\sin\frac x2)$ and thus setting $$ 2\sin\frac{x(t)}2 = 2\sin\frac 12\cdot \sin(u(t)),\\ \dot x(t)= 2\sin\frac 12\cdot \cos(u(t)), $$ $u(0)=\frac\pi2$, and then comparing the derivative of the first equation with the second equation, $$ \cos\frac{x(t)}2\cdot\dot x(t) = 2\sin\frac 12\cdot\cos(u(t))\dot u(t) $$ so that for non-constant solutions (so that $\cos(u(t))\ne 0$) $$ \dot u(t) = \cos(\frac{x(t)}2)=\sqrt{1-\sin^2\frac 12\cdot \sin^2(u(t))} $$ For the plot you have then to reverse the parametrization to $$ x(t)=2\arcsin\left(\sin\frac 12\cdot \sin(u(t))\right) $$

9
On

$$ \frac{d^2x}{dt^2} + \sin(x) = 0 \tag 1$$ $$ \frac{d^2x}{dt^2}\frac{dx}{dt} + \sin(x)\frac{dx}{dt} = 0\tag 2$$ When you multiply Eq.$(1)$ by $\frac{dx}{dt}$ in order to get Eq.$(2)$, you introduces additional solutions which are solutions of $(2)$ but not solutions of $(1)$. The additional solutions are the solutions of $$\frac{dx}{dt}=0 \quad\implies\quad x(t)=c\tag 3$$ The boundary conditions $x(0)=1$ and $x'(0)=0$ give $c=1$. So an additional solution was introduced : $$x(t)=1$$ Thus is is correct that, according to the boundary conditions, the drawing of the solution of $(2)$ is the same as the drawing of the solution of $(1)$ plus the straight line $x=1$.

$$ $$

IN ADDITION after the comments :

The analytic solution of $\quad\frac{dx}{dt}=\sqrt{2(\cos(x)-\cos(1))}\quad$ is $$\begin{cases} x(t)=1\\ x(t)=2\text{ am}\left(\sqrt{\frac{1-\cos(1)}{2}}(c-t) \:\bigg|\: \csc^2(\frac12) \right) \end{cases}$$ $\text{am}$ denotes the Jacobi amplitude function.

With the condition $x(0)=1\quad\implies\quad c=\csc(\frac12)F\left(\frac12 \:\big|\:\csc^2(\frac12) \right)\simeq 1.67499$

$F$ denotes the elliptic integral of the first kind.

enter image description here

2
On

Apparently the preceding answers don't satisfy you because they are more mathematical answers than about software and/or about how to handle software (Matlab).

I am not matlab expert. So my answer cannot be considered as definitive and might be discussed.

You coded the equation $$x'(t)=\frac{dx}{dt}=\sqrt{2\left(\cos(x)-\cos(1)\right)} \tag 1$$ With the initial specification $x(0)=1$.

From the behaviour of software that you report, one can suppose that $x'_0=x'(0)=\sqrt{2\left(\cos(1)-\cos(1)\right)} =0$ was first computed.

Then applying a time increment $\delta t$ at time $t=0$ the nest time $t_1=0+\delta t= \delta t$ and the next point $(x_1,t_1)$ was computed : $x_1=x_0+x'_0\delta t=1+0(\delta t)=1$.

And so on for successive increments, $x'$ remains nul and $x$ remains equal to $1$. So the software give a correct solution $$x(t)=1.$$ As already pointed out, this is not the unique solution of the ODE Eq.$(1)$ with condition $x(0)=1$.

Apparently the software doesn't give spontaneously another solution. As a consequence one have to code the software in order to get another solution since we know that several solutions exist.

Consider the first step of the numerical calculus : $\begin{cases} t_1=\delta t \\ x_1=1+\delta x_1 \end{cases}$

In Eq.$(1)$ necessarily $\quad \cos(x_1)-\cos(1)> 0$ thus $x_1<1$ and $\delta x_1<0$.

$\frac{\delta x_1}{\delta t}<0 \quad$ is not consistent with Eq.$(1)$ where the square root is positive.

Thus Eq.$(1)$ is not correct. The correct equation according to the initial condition $x(0)=1$ is $$x'(t)=\frac{dx}{dt}=-\sqrt{2\left(\cos(x)-\cos(1)\right)} \tag 2$$ This is a correct derivation of $\frac12(\frac{dx}{dt})^2=\cos(x)-\cos(1)$.

$\underline{\text{So, first you have to correct the code}}$ to apply Eq.$(2)$ instead of $(1)$.

But how to make the algorithm start a non-constant incremental process ? Preliminary inspection shows that :

$\delta x=-\sqrt{2\left(\cos(1+\delta x)-\cos(1)\right)}\:\delta t=-\sqrt{2\left(\cos(1)\sin(\delta x))-\sin(1)\sin(\delta x)-\cos(1)\right)}\:\delta t\simeq -\sqrt{-2\sin(1)\delta x}\:\delta t$

After simplification :$\quad \delta x=-2\sin(1)(\delta t)^2$

The first increment $\delta t$ gives : $$\begin{cases} t_1=\delta t \\ x_1=1-2\sin(1)(\delta t)^2 \end{cases}$$

$\underline{\text{Second, you have to correct the code }}$ for the above starting point instead of $t=0\:,\:x=1$.

You mentioned in comment, citation : "I already did that with 0.9 and 0.8. Instead on line, I am getting something else but not sinusoidal curve". But this was not exactly the right way.

First, change the sign of the square root. Second, try with a smaller $\delta t$, for example : $$\begin{cases} t_1=0.001 \\ x_1=1-2\sin(1)(0.001)^2 \end{cases}$$