Introduction to MATLAB Exercise 4

Brocard's Problem. It is conjectured by the Paul Erdos that there are only three pairs of \( (m,n)\) to the following equation \[ n! + 1= m^2 \] for positive integer \(m\) and \(n\). Find this three solutions using a loop for \(n\) taking values from \(1\) to \(17\) (MATLAB can not handle larger integers like \(18!\), and you may have to use symbolic integer that will be covered next week). Think about how to make your implementation less error-prone, when checking the condition \(n! + 1= m^2\).
Wilson prime. The famous Wilson's theorem states that a number \(n\) is prime, if and only if \[ (n-1)!\ \equiv\ -1 \ \ (\mbox{mod } n). \] That is \( (n-1)!+1\) is divisible by \(n\). You can verify the first few integers like \(n=2,3,4,5\), by using (change the value of \(n\))
```
n = 2; mod(factorial(n-1)+1,n) % you get zero, if n is a prime; otherwise, you get a non-zero number
```
There are a few special prime numbers \(p\) called Wilson prime, besides the condition \( (p-1)!\ \equiv\ -1 \ \ (\mbox{mod } p)\), which satisfy the more restricted relation \( (p-1)!\ \equiv\ -1 \ \ (\mbox{mod } p^2)\). The first is \(p=5\): frome \( (5-1)! = 4! =24\), we get \( (5-1)! \equiv -1 (\mbox{mod } 25)\), instead of just \( (5-1)! \equiv -1 (\mbox{mod } 5)\). Complete the following script, to find all such Wilson primes that are less than 1000 (there is only one between 10 and 100, and then only another one between 100 and 1000).
```
N = 1000;
pl = primes(N);  % get the list of all primes less than or equal to N
for p=pl  % loop over all primes up to N
    % check whether p is a Wilson prime
end
```
Because the built-in function factorial can only be used for numbers less than 19 (forget about the much slower alternative using symbolic numbers), you have to apply the modulus frequent enough.
Solve the ODE \begin{equation} y''=\frac{y'^2}{y}-y^2. \end{equation} symbolically using dsolve, whose general solution is \begin{equation} y(x) = \frac{1-\tanh^2\frac{x+c_2}{2c_1}}{2c_1^2}. \end{equation} You can use diff(y,x,2) for the second derivative of y with respective to x, and similarly for higher order derivatives. Notice that the constants in the general solution could be different from the above expression (but still equivalent).
Solve the ODE \begin{equation} y'(x) = \frac{3y(x)}{x}+\frac{x^5}{2y(x)+x^3}, \end{equation} symbolically with the initial condition \(y(1)=0\), which should lead to the exact solution \begin{equation} y(x) = \frac{\sqrt{1+4\ln x}-1}{2}x^3. \end{equation} Now try to find the solution on the interval \([1,2]\) using ode45, and compare the error between the numerical solution and the exact solution (by evaluating at the vector of time from the numerical solver).
(Chaotic path driven by periodic functions) Plot the trajectory of the following differential equation on the phase place, that is the \( (x, \dot{x})\) plane: \begin{equation} \ddot{x} + 3\dot{x}-4|\dot{x}|+x=\sin t,\qquad x(0)=1,\ \dot{x}(0)=-0.2, \ t \in [0, 500]. \end{equation} You may define options by adding one line options = odeset('refine',8); in your code, and then ode45( ..., options) to get enough points/resolution for plotting. This system is chaotic, because the solution could be completely different in the long run, if the initial condition is changed slightly. Compare the (numerical) solution at the end of the time interval \(t=500\) with the initial condition \(x(0)=1,\dot{x}(0)=-0.2\) and the initial condition \(x(0)=1+10^{-10},\dot{x}(0)=-0.2\).
(Chaotic path driven by periodic functions) Show the trajectory of the following differential equation on the phase place, that is the \( (x, \dot{x})\) plane: \begin{equation} \ddot{x} = (x-\dot{x})\sin t - (x+1)\dot{x},\qquad x(0)=1,\ \dot{x}(0)=0.5, \ t \in [0, 500]. \end{equation} You may define options by adding one line options = odeset('refine',8); in your code, and then ode45( ..., options) to get enough points/resolution for plotting. You can also change the initial condition slightly, to see how the final solution at \(t=500\) changes drastically.
The non-linear ordinary differential equation \begin{equation} y' = -y^3 + \sin t \end{equation} has a periodic solution with period \(2\pi\) (that is \(y(t+2\pi)=y(t)\) for any \(t\geq 0\)) for some intial condition \(y_0=y(0)\). From the existence and uniqueness of solutions of this simple ODE, you only need to make sure \(y(2\pi)=y(0)\), to get such a periodic solution. Use this condition to find the right initial condition \(y_0\). Hint: Solve the ODE numerically on the interval \([0,2\pi]\), and find \(y_0\) such that it is a root of \(y(2\pi)-y_0=0\).
```
      
%% Periodic solution of the non-linear ODE y' = -y^3 + sin (t)
% save this file as per_nlode.m
function per_nlode
y0 = fzero(@(y0)sol_diff(y0),0); 

function yd = sol_diff(y0)
% return the difference y(2*pi)-y(0) of solutions to y' = -y^3 + \sin t, with input y(0)=y0
```
Once you find the right initial condition \(y(0)\), solve it on the interval \([0,20\pi]\) and plot it (which should look periodic with period \(2\pi\)).