[Math] How to solve the ODE using Newton’s method

Question

Let consider on $[0,20]$, the ODE : $y'(t)=y^2(t)-y^3(t)$ and $y(0)=1$.

I start to use a backward Euler method using : $y_{n+1}=y_n+h(y_{n+1}^2-y_{n+1}^3)$ with $h>0$ given.

I build $F:y_{n+1}\mapsto y_{n+1}-y_n-h(y_{n+1}^2-y_{n+1}^3)$ such that $F(y_{n+1})=0$.

If $F'(y_{n+1})\neq 0$, I can use the Newton's method to determine $y_{n+1}$.

Do I have to consider $y_{n+1}-\frac{F(y_{n+1})}{F'(y_{n+1})}=y_{n+1}-\frac{y_{n+1}-y_n-h(y_{n+1}^2-y_{n+1}^3)}{1-2hy_{n+1}+3hy_{n+1}^2}=\frac{y_n-hy_{n+1}^2+2hy_{n+1}^3}{1-2hy_{n+1}+3hy_{n+1}^2}$ ? Then how to continue ?

Thanks in advance !

Answer 1

obviously your ODE has the stationary solution $y(t)=1$.

Anyway,your ODE is clearly non-linear, so, as you stated, you have to set to zero the following function:

$F(x)=x-y_n - hf(x)$, where $f$ is defined as $ f=y^2 - y^3$.

So, $F'(x)= 1 - h\frac{df}{dx}(x)$ (i used this notation because $f$ could be a function of $t$)

In order to do this, you have to use Newton's method:

given $x_1=y_n$ (the current value of the solution is the initial guess for Newton's iteration), do $x_{k+1}=x_k - \frac{F(x_k)}{F'(x_k)}$ until the difference $|x_{k+1} - x_k|$ or the norm of the 'residue' is less than a given tolerance (or combination of absolute and relative tolerances)

This can be coded as follow. In this code I used Newton's method for non-linear systems.

clear all
close all
f=@(y) y^2 - y^3;
df=@(y) 2*y - 3*y^2; %derivative of f (in a systems of ODE, it will be 
%the jacobian matrix of f)

ts=1000; %time steps
k=20/ts;
y0=1; %initial condition

%function to set to zero:
%y will be the solution at time n+1   
%yn is the current y
F=@(y,yn) y-yn-k*f(y);
J=@(y) 1 - k*df(y);%jacobian matrix for F. In this case is just a number

y=NaN(1,ts+1);%initialize the solution to NaN
y(1,1)=y0;
tol=k^2; %using backward Euler

for n=1:ts
    y(1,n+1)=y(1,n); %initial guess for Newton's method

    res=-J(y(1,n+1))\F(y(1,n+1),y(1,n)); %Deltax
    while (norm(res,inf)>tol)
        y(1,n+1)+=res;
        res=-J(y(1,n+1))\F(y(1,n+1),y(1,n));
    end
    y(1,n+1)+=res;
 end

t=linspace(0,20,ts+1);
plot(t,y,'k--')
xlabel('t')
ylabel('y')

And the result agrees with the previous calculation, as shown in the graphic.

Please notice that I set the 'tol' parameter proportional to $h^2$: this is due to the fact that the Backward Euler's error is $\mathcal{O}(k^2)$