numerical methods
How can I show that the Euler method fails to approximate the exact solution
$$y(x)=(2x/3)^{3/2}$$
to the IVP
$$y'=y^\frac{1}{3}$$
$$y(0)=0$$
Here we have $f(t,y)= y^\frac{1}{3}$, $y_0=0$ and so $f(t_0,y_0)=f(0,0)=0$ and
$$y_{n+1}=y_n +h f(t_n,y_n)$$
Thus
$$y_1=0 \\
y_2=0 \\
\vdots\\
y_n=0$$
So, I can't understand why it fails. Could you help me?
When you implement the implicit Euler method, the equation to be solved at every time step is not solved in an exact way, but rather by some numerical method, typically the fixed point method. You can show that for nice enough $f$ you can get convergence for small $h$, namely $h <\frac 1L$ (where $L$ is the Lipschitz constant for $f$). The method would read something like
$$ \begin{cases} y_0\\ \begin{cases} y_{n+1}^{(0)} = y_n\\ y_{n+1}^{(k+1)} = y_n + h f(t_{n+1},y_{n+1}^{(k)}), \quad k \ge 0 \end{cases} \end{cases}. $$
This way (with this choice of initial approximations), the method keeps selecting solutions "close" to the solution obtained in the previous time step.
$y_{n+1}$ is obtained as the solution of $y_{n+1} = y_n + h f(t_n+h, y_{n+1})$, which is an equation of the form $y_{n+1} = \phi(y_{n+1})$, i.e. $y_{n+1}$ is a fixed point of $\phi$. The natural way to solve this equation is to use the fixed point method, that starts with an initial approximation of $y_{n+1}$, that I denoted as $y_{n+1}^{(0)}$ and builds a sequence $y_{n+1}^{(k+1)} = \phi(y_{n+1}^{(k)})$ with the property $$ \lim_{k\to \infty}y_{n+1}^{(k)} = y_{n+1}. $$
You can avoid many index errors if you construct the occurring array variables via array operations. Here one can compress the code to
N=100;
ta=1;
tb=exp(1);
h=(tb-ta)/N;
x=linspace(ta,tb,N+1);
% y''+4/x*y'-2/x^2*y = -2/x^2*ln(x)
yeks = x.^(-3/2+sqrt(17)/2)+x.^(-3/2-sqrt(17)/2)+3/2+log(x);
n=length(x);
p = 4 ./ x;
q = -2 ./ x.^2;
a = 1/h^2-p ./ (2*h);
b = q-2/h^2;
c = 1/h^2+p ./ (2*h);
%Coefficient matrix for the equations for x(2) .. x(N)
A=spdiags([ a(3:N+1); b(2:N); c(1:N-1) ]', [-1 0 1], N-1, N-1);
%RHS matrix
B=-2 ./ x.^2.*log(x);
B = B(2:N)';
B(1) -= a(2)*yeks(1);
B(N-1) -= c(N)*yeks(N+1);
ynum = [ yeks(1); A\B; yeks(N+1) ];
One can use several diag commands to construct the tridiagonal matrix from its sub-diagonals, or use the sparse banded matrix constructor spdiags as above. This gives an error level of $10^{-5}$ and the plot
Best Answer
Note that the general solution, $$ y(x) = \left(\frac{2x+C}{3}\right)^{3/2}$$ Can be achieved by integration: $$\int \frac{dy }{y^{1/3}}=\int dx$$ This assumes that $y\neq 0$, which is not necessarily true. You therefore have two solutions at $x=0$, and you can't force the Euler method to "choose" the right one.
The reason you don't have a unique solution at $x=0$ is that $d(y^{1/3})/dx$ does not exist on any open interval containing $x=0$. Choosing an IV point with $x> 0$ will guarantee local uniqueness.