Newton's first differential equation

One of the earliest known differential equations, which Newton solved with series expansion in 1671 is: $$ \begin{align} y'(x) =1-3x+y+x^2+xy,\ y(0)=0 \label{eq:newton} \end{align} $$

Newton gave the following solution: $$ \begin{align} y(x) \approx x-x^2+\frac{x^3}{3}-\frac{x^4}{6}+ \frac{x^5}{30}-\frac{x^6}{45} \label{eq:newton_approx} \end{align} $$

Today it is possible to give the solution on closed form with known functions as follows, $$ \begin{align} y(x)=&3\sqrt{2\pi e}\cdot \exp\left[x\left(1+\frac{x}{2}\right)\right]\cdot \left[\text{erf}\left(\frac{\sqrt{2}}{2}(1+x)\right)-\text{erf}\left(\frac{\sqrt{2}}{2}\right)\right] + 4\cdot\left[1-\exp[x\left(1+\frac{x}{2}\right)\right]-x \label{eq:newton_analytical} \end{align} $$

Note the combination \( \sqrt{2\pi e} \).

Problem 1: Euler's method

Solve Eq. \eqref{eq:newton} using Euler's method. Plot and compare with Newton's solution Eq. \eqref{eq:newton_approx} and the analytical solution Eq. \eqref{eq:newton_analytical}. Plot between \( x=0 \) and \( x=1.5 \)

Hint 1.

Hint 2.

Mathematical pendulum with friction

Figure 2: Mathematical pendulum.

Figure 2 shows a mathematical pendulum. By modeling the friction forces acting on the pendulum as \( \mu \frac{d \theta}{dt} \), the motion can be described by the following ODE: $$ \begin{align} \frac{d ^2\theta}{dt^2} + \frac{\mu}{m} \frac{d \theta}{dt} + \frac{g}{l} \sin\theta &=0, \label{eq:1} \end{align} $$

where \( \mu \) is a friction coefficient, \( l \) is the length og the rod and \( m \) is the mass of the pendulum with initial conditions $$ \begin{align} &\theta(0) =\theta_0 \label{_auto1}\\ &\frac{d \theta}{dt}(0)=\dot\theta_0. \label{_auto2} \end{align} $$

Assume that the pendulum wire is a massless rod, such that \( -\pi \leq \theta_0 \leq \pi \). Unless stated otherwise, use \( \mu=1 \), \( m=1 \), \( g=9.81 \), \( l=1 \), \( \theta_0=85^o \) and \( \dot\theta_0=0 \).

Problem 2: Euler's method

Write a Python program that solves the ODE in \eqref{eq:1} with the specified initial conditions using Euler's method (remember to use radians). Experiment with different time steps \( \Delta t \), and friction coefficients \( \mu \). How do they affect the stability of the solution? Carry out the computation for a few cycles, and plot the amplitude as a functions of \( t \).

Problem 3: Heun's method

a) Solve Problem 2 using Heun's method.

b) Write a function that detects the amplitude peaks, and use this to calculate the period of the pendulum. Is it constant?

Problem 4: (Optional) General second order scheme

a) Solve Problem 2 using Midpoint method and Ralston's method. See: "See http://folk.ntnu.no/fredf/teaching/tkt4140/exercises/pythonExercises/pexercise2/TKT4140%20Lec%204%20RKmethods%20(1).pdf

c) create an animation similar to the one above.