Accuracy of Explicit Euler method (finite difference) decreases as Δx decreases, shouldn't it increase?

Question

The price of a commodity can be described by the Schwartz mean reverting SDE $$dS = \alpha(\mu-\log S)Sdt + \sigma S dW, \qquad \begin{array}.W = \text{ Standard Brownian motion} \\ \alpha = \text{ strength of mean reversion}\end{array}$$

From it is possible to derive the PDE for the price of the forward contract having the commodity as underlying asset $$\tag1\frac{\partial F}{\partial t} + \alpha\Big(\mu-\frac{\mu-r}\alpha -\log S\Big)S\frac{\partial F}{\partial S}+\frac12\sigma^2S^2\frac{\partial^2F}{\partial S^2} = 0$$

whose analytical solution is

$$F(S,\tau)=\exp\bigg(e^{-\alpha\tau}\log S +\Big(\mu-\frac{\sigma^2}{2\alpha}-\frac{\mu-r}{\alpha}\Big)(1-e^{-\alpha\tau})+\frac{\sigma^2}{4\alpha}(1-e^{-2\alpha\tau})\bigg)$$

where $\tau=T-t$ is the time to expiry ($T$ is the time of delivery/expiry).

Using Euler explicit method, i.e. forward difference on $\dfrac{\partial F}{\partial t}$ and central difference on $\dfrac{\partial F}{\partial S}$ and $\dfrac{\partial^2F}{\partial S^2}$, we can discretize eq (1) as $$F^{n+1}_i = a F^n_{i-1} + b F^n_i + c F^n_{i+1}$$ where $a = \dfrac{S\Delta t}{2\Delta S}\bigg(\alpha\mu-(\mu-r)-\alpha\log(S)-\dfrac{\sigma^2S}{\Delta S}\bigg)$

$b = \bigg(1-\sigma^2S^2\dfrac{\Delta t}{\Delta S^2}\bigg)$ and $c = \dfrac{S\Delta t}{2\Delta S}\bigg(-\alpha\mu+(\mu-r)+\alpha\log(S)-\dfrac{\sigma^2S}{\Delta S}\bigg)$.

To run Explicit Euler we have then to choose the number $N$ of time steps, which also set $\Delta t$ since $\Delta t = T/N$, and the size of $\Delta S$. Since the finite difference scheme divides the cartesian plane (time is the X-axis, and spot price is the Y-axis) in a grid, if we take more time steps and/or smaller $\Delta S$ the grid will be more dense and the accuracy of the approximation should increase.

However, the code I wrote using the equations above doesn't work in this way, in particular to have big accuracy I have to use a large $\Delta S$, and the accuracy decreses when using small values of $\Delta S$ to point that by using $\Delta S=0.1$ the relative error explodes to $10^{165}$ as you can see in the image below (dS stands for $\Delta S$).

Here is the matlab code, I think there is an error somewhere, moreover I'm not sure about the boundary conditions F(1) and F(end)

%% Data and parameters
spot_prices = [ 22.93 15.45 12.61 12.84 15.38 13.43 11.58 15.10 14.87 14.90 15.22 16.11 18.65 17.75 18.30 18.68 19.44 20.07 21.34 20.31 19.53 19.86 18.85 17.27 17.13 16.80 16.20 17.86 17.42 16.53 15.50 15.52 14.54 13.77 14.14 16.38 18.02 17.94 19.48 21.07 20.12 20.05 19.78 18.58 19.59 20.10 19.86 21.10 22.86 22.11 20.39 18.43 18.20 16.70 18.45 27.31 33.51 36.04 32.33 27.28 25.23 20.48 19.90 20.83 21.23 20.19 21.40 21.69 21.89 23.23 22.46 19.50 18.79 19.01 18.92 20.23 20.98 22.38 21.78 21.34 21.88 21.69 20.34 19.41 19.03 20.09 20.32 20.25 19.95 19.09 17.89 18.01 17.50 18.15 16.61 14.51 15.03 14.78 14.68 16.42 17.89 19.06 19.65 18.38 17.45 17.72 18.07 17.16 18.04 18.57 18.54 19.90 19.74 18.45 17.33 18.02 18.23 17.43 17.99 19.03 18.85 19.09 21.33 23.50 21.17 20.42 21.30 21.90 23.97 24.88 23.71 25.23 25.13 22.18 20.97 19.70 20.82 19.26 19.66 19.95 19.80 21.33 20.19 18.33 16.72 16.06 15.12 15.35 14.91 13.72 14.17 13.47 15.03 14.46 13.00 11.35 12.51 12.01 14.68 17.31 17.72 17.92 20.10 21.28 23.80 22.69 25.00 26.10 27.26 29.37 29.84 25.72 28.79 31.82 29.70 31.26 33.88 33.11 34.42 28.44 29.59 29.61 27.24 27.49 28.63 27.60 26.42 27.37 26.20 22.17 19.64 19.39 19.71 20.72 24.53 26.18 27.04 25.52 26.97 28.39 ];
S = spot_prices; % real data

r = 0.1;    % yearly instantaneous interest rate
T = 1/2;   % expiry time

alpha = 0.069217; %
sigma = 0.087598; % values estimated from data
mu = 3.058244;    %

%% Exact solution
t = linspace(0,T,numel(S));
tau = T-t; % needed in order to get the analytical solution (can be seen as changing the direction of time)
F = exp( exp(-alpha*tau).*log(S) + (mu-sigma^2/2/alpha-(mu-r)/alpha)*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) ); % analytical solution
F(1) = 0; % I think since there is no cost in entering a forward contract
plot(t,S)
hold on
plot(t,F,'g')
Exact_solution = F;

%% Explicit Euler approximation of the solution
S1 = S(2:end-1);  % all but endpoints
N = 3000; % number of time steps
dt = T/N; % delta t
dS = 1e1; % delta S, by decreasing dt and/or dS the approximation should improve
for m = 1:N
    F(2:end-1) = S1*dt/2/dS.*( alpha*mu-(mu-r)-alpha*log(S1)-sigma^2*S1/dS).*F(1:end-2) ...
               +                                  (1+sigma^2*S1.^2*dt/dS^2).*F(2:end-1) ...
               + S1*dt/2/dS.*(-alpha*mu+(mu-r)+alpha*log(S1)-sigma^2*S1/dS).*F(3:end);
    F(1) = 0; % correct?
    F(end) = S(end); % correct?
end
plot(t,F,'r.')
legend('Spot prices','Forward prices from exact solution','Forward prices from Explicit Euler')
title("dS = " + dS + ", relative error = " + norm( F-Exact_solution,2 ) / norm( Exact_solution,2 ))
xlabel('time')
ylabel('price')

Answer 1

I haven't inspect everything in detail but I think I have a useful remark. Since you know the analytical solution, you can compute from there the boundary conditions.

For $t = T$, you have $\tau = 0$. You get $F(t = T; \tau = 0) = S(T)$ when replacing in the analytical solution. Then, F(end) = S(end), as you stated in your code.

On the other hand, when $t = 0$, $\tau = T$. If you replace that in the analytical solution, you get the same expression that you wrote but replacing $S$ with $S(0)$ and $\tau$ with $T$ so, in your code, F(1) = 0 do not match with the analytical solution. Since I also believe that F(1) should be zero, there might be something that we are missing in the analytical solution...

Please, forgive me if I am mistaken, but I thought it would be a useful comment.

PS: comments below address the underlying problem.