First-order ODE scheme implementation giving less than first-order convergence?

Question

I am solving the initial value problem $$ \frac{d}{dt} (E C_g) = -\delta, \quad E(0) = E_0, $$ for $E$, where $E$ and $C_g$ are functions of $t$, $C_g$ is completely known, and $\delta$ is a function of $E$. Using the following first order approximations for derivatives, \begin{align*} \frac{d}{dt} (E C_g) = E' C_g + E C_g', \\ E' \approx \frac{E_i - E_{i-1}}{\Delta t}, \quad C_g' \approx \frac{C_{gi} - C_{gi-1}}{\Delta t}, \end{align*} we end up with the following finite difference scheme to step forward in time \begin{align*} E_i = \frac{-\delta_{i-1} \Delta t + 2E_{i-1} C_{gi-1} - E_{i-1}C_{gi}}{C_{gi-1}}. \end{align*}

When I implement this in MATLAB, I test my code by solving this IVP with a specific $E$, which leads to a specific right-hand side, and then can compute my error by comparing my computed solution to my true solution using

 error = norm(E_computed - E_true).

I can then refine the grid spacing and save the spacing values in "spacingArray" and the error values in "errorArray", then I can estimate the order of convergence using

 order = (log(errorArray(3)) - log(errorArray(6)))/(log(spacingArray(3)) - log(spacingArray(6))),

but when I do this, I get order $\approx 0.5$, instead of order $\approx 1$ (this can also be seen by looking at a log-log plot of spacing against error). Does this mean that there is something wrong with the way that I have implemented the finite difference scheme in MATLAB, or could there be something else wrong?

Since my problem is a little more nuanced than I originally explained, here is the exact MATLAB code I run:

%% Wave-energy model verification code for d/dx (E*C_g) = -delta
% To verify the model, we apriori choose H, h, and sigma;
% since d/dx (E*C_g) + delta = err, where err is some error.
% We obtain an analytical formula for err by determining
% d/dx (E*C_g) + delta (by using Mathematica or pencil/paper).
% Then, we solve the differential equation d/dx (E*C_g) = -delta + err,
% where the delta and err depend on our apriori chosen values of H, h,
% and sigma. 
% 7/2/2017 9:38 AM

%% Loop indices
start = 3; jump = 8;

%% Initial arrays for testing
errorArray = zeros(jump, 1);
dxArray = zeros(jump, 1);
numPointArray = zeros(jump, 1);

for kk = start:(start + jump)

    %% Set up the "mesh"
    domain = [0,10];
    numPoint = 2^kk;
    numPointArray(kk - start + 1) = numPoint;
    x = linspace(domain(1), domain(2), numPoint);
    dx = (domain(2) - domain(1))/(numPoint - 1);
    dxArray(kk - start + 1) = dx;

    %% Physical parameters
    rho = 1000; % water density
    g = 9.8;    % gravity 
    gam = .78;  % ?
    rms = .707; % root-mean square constant
    B = 1;      % energy dissipation parameter


    %% Apriori choices for h, H, and sigma:
    h = @(x) 50 - x;
    h = h(x);
    h_prime = @(x) -1*ones(size(x));
    h_prime = h_prime(x);
    sigma = @(x) 2*ones(size(x));
    sigma = sigma(x); f = sigma/(2*pi);
    Ha = @(x) cos(x) + 5; % Ha = "analytical height"
    Ha = Ha(x);
    Ha_prime = @(x) -sin(x);
    Ha_prime = Ha_prime(x);


    %% Solving for wavenumber via the "linear" dispersion relation
    dispersionRelation = @(k) sigma.^2 - g.*k.*tanh(k.*h);
    k = ones(size(x)); % random initial guess for iterative solver
    k = fsolve(dispersionRelation, k); 

    %% Off-shore boundary conditions
    H = zeros(size(x)); E = zeros(size(x)); H_rms = zeros(size(x));
    H(1) = 6; % offshore wave height determines H_rms(1) and E(1)
    H_rms(1) = rms*H(1);        
    E(1) = (1/8)*rho*g*H_rms(1).^2;
    R = zeros(size(x)); delta = zeros(size(x));
    H_b = gam*h;
    R(1) = H_b(1)/H_rms(1);

    %% Group Celerity
    c = sigma./k; 
    C_g = (c/2).*(1 + (2*k.*h)/(sinh(2*k.*h)));

    %% Energy flux, F = E*C_g
    F = zeros(size(x));
    F(1) = E(1)*C_g(1); % boundary condition

    %% -------------------------------------------------------------------
    %% The analytical formulae (plural) for d/dx (E*C_g) from Mathematica
    % using shallow water approximation for wavenumber:
    theta1 = 2*sigma.*sqrt(h)/sqrt(g);
    shallowFluxDerivative = (1/16)*rho*g^(3/2)*(rms)^2*(Ha.^2.*(2.*sigma.*csch(theta1).* ...
        h_prime/sqrt(g) + h_prime./(2*sqrt(h)) - 2*sigma.^2.*coth(theta1).* ...
        csch(theta1).*sqrt(h).*h_prime/g) + 2*(sqrt(h) + 2*sigma.*csch(theta1).*h/sqrt(g)).* ...
        Ha.*Ha_prime);
    % using the deep water approximation for wavenumber:
    theta2 = 2*h.*sigma.^2/g;
    deepFluxDerivative = (1/16)*rho*g^2./sigma.*(rms)^2.*(Ha.^2.* (2*sigma.^2.*csch(theta2).*h_prime/g - ...
        4*sigma.^4.*coth(theta2).*csch(theta2).*h.*h_prime/(g^2) ) + 2*(1 + 2*sigma.^2.*csch(theta2).*h/g).* ...
        Ha.*Ha_prime);

    %% Analytical formula for delta
    Ha_rms = rms*Ha;
    Ra = H_b./Ha_rms;
    deltaAnalytic = (1./(4.*h)).*B.*rho.*g.*(sigma./(2.*pi)).*Ha_rms.^3.*( (Ra.^3 + 1.5.*Ra).*exp(-Ra.^2) + 0.75.*sqrt(pi).*(1 - erf(Ra)));

    %% Analytical "err" formula, d/dx(flux) + delta = err
    shallowError = shallowFluxDerivative + deltaAnalytic;
    deepError = deepFluxDerivative + deltaAnalytic;
    %% -------------------------------------------------------------------

    %% Solve the ODE via finite differences
    for ii = 2:length(x)
        F(ii) = -delta(ii - 1)*dx + F(ii - 1); % finite difference formula
                                               % to compute new flux

        % determine which regime we are in, then add the approriate error
        if h < 0.05*2*pi./k(ii) % shallow
            F(ii) = F(ii) + dx*shallowError(ii);
        elseif h > 0.5*2*pi./k(ii) % deep
            F(ii) = F(ii) + dx*deepError(ii);
        else % default to deep?...
            F(ii) = F(ii) + dx*deepError(ii);
        end

        E(ii) = F(ii)./C_g(ii); % compute new energy from flux
        H(ii) = sqrt(8./(rho*g)*E(ii))/rms; % compute height from energy
        H_rms(ii) = rms*H(ii); 
        delta(ii) = 1/(4*h(ii))*B*rho*g*f(ii)*H_rms(ii)^3*(R(ii)^3 + (3/2)*R(ii))*exp(-R(ii)^2) + ...
            (3/4)*sqrt(pi)*(1 - erf(R(ii))); % compute new delta
    end

    %% post processing
    error = norm(H - Ha); % difference between computed and true solution
    errorArray(kk - start + 1) = error;

end

%% create a log-log plot
figure(), loglog(dxArray, errorArray), title('log-log plot of error against spacing')
axis equal, xlabel('log(dx)'), ylabel('log(norm(abs(H - H_a)))')

%% determine the slope of the log-log graph 
% this should match with the theoretical order of convergence
order = (log(errorArray(jump)) - log(errorArray(1)))/(log(dxArray(jump)) - log(dxArray(1)))

For some reason, I can't save and upload the plot, so here is the outputted order of convergence: 0.512087039390502.

One thing that is interesting is that it looks like the order of convergence is a little bit better for the first refinement, and in fact I can compute

 order = (log(errorArray(2)) - log(errorArray(1)))/(log(dxArray(2)) - log(dxArray(1)))

which produces $\text{order} \approx 0.629469265412451$, which is at least a little closer...

Answer 1

Let me modify slighly your ODE to see the essence of the problem. Let $$\left\{\begin{array}{l}y'=-y \\ y(0) = 1 \end{array} \right. \qquad t\in[0,T]\tag{*}$$ the ODE to be solved. Everyone knows that the exact solution of $(*)$ is: $$y=\exp{(-t)} \tag{*a}$$

If we discretise $(*)$ one would obtain: $$\frac{y^{n+1}-y^n}{\Delta t}=-y^{n+1} \qquad y^0=1\tag{*b}$$ in which have followed the unconditionally stable scheme being the RHS of $(*)$ implicitly treated.

In fact equation $(*b)$ solves with second order of accuracy the following problem: $$\left\{\begin{array}{l}y'=-y-\Delta t\,y'' \\ y(0) = 1 \end{array} \right.$$

The exact numerical solution of $(*b)$ is: $$y^n=(1+\Delta t)^{-n}$$ At time $t_i$, with $i\in[0,N]$ where $N$ is large wrt. $T$, the error would be $$\epsilon_i=|y(t_i)-y^{i}|=|\exp{(-t_i)}-(1+\Delta t)^{-i}|=|\exp{(-iT/N)}-(1+T/N)^{-i}|\approx \frac{iT}{2N}\frac{T}{N}=\frac{t_i}{2}\Delta t=\mathcal{O}(\Delta t) $$

Therefore, we can conclude from the latter equation that the error can be bounded by a constant that depends only on the domain in which $(*)$ is integrated, $i.e.$ $$\epsilon\leq C(T)\Delta t=T\,\Delta t \tag{*c}$$

This analysis can be performed with any ODE leading to the same order in $\Delta t$ according to the chosen scheme.

The conclusion $(*c)$ can be determined numerically if we notice that the error in the $i$th time is bounded by the same power $p$ of $\Delta t$ with the same proportionality constant, $i.e$ for two different meshes $\Delta t$ and $\Delta t'$ we will obtain errors $\epsilon$ and $\epsilon'$ respectively such that $$\frac{\epsilon'_i}{\epsilon_i}=\left(\frac{\Delta t'}{\Delta t}\right)^{p} \tag{*d}$$

To use $(*d)$ choose, for example the last node (which always independent on the used mesh and it coincides with $t_i=T$) and calculate $p$ with $$p=\frac{\log{\left[\epsilon'(T)/\epsilon(T)\right]}}{\log{\left[\Delta t'/\Delta t\right]}}$$

This must be always around $p=1$ and is always dependent on the discretisation method and the regularity properties of the $RHS$ (derivatives) of $(*)$ that can negatively interfere in the order of the method.

What is more I leave you a code that summarises what I previously said:

T = 1;

deltat1 = 1e-3; %Multiple of 10
deltat2 = 1e-4; %Multiple of 10

%Exact solution evaluated at $T$
yTexact = exp(-T);

%Numerical solution for the 1st mesh
N = T/deltat1;
y1 = (1+deltat1).^(-[0:N]);
y1T = y1(end); %Numerical solution at t=T
epsilon1 = abs(yTexact-y1T); %error for the last node

%Numerical solution for the 2nd mesh
N = T/deltat2;
y2 = (1+deltat2).^(-[0:N]);
y2T = y2(end); %Numerical solution at t=T
epsilon2 = abs(yTexact-y2T); %error for the last node

%Compute the method order:
p = log(epsilon2/epsilon1)/log(deltat2/deltat1)

The result is $p=0.998$

EDIT If you compute the global error, $i.e.$ the norm (for $N$ large): $$\sqrt{\sum{\epsilon_i^2}}\approx \frac{T}{2N}\frac{T}{N}\sqrt{\sum_{i=1}^{N}{i^2}}\approx \frac{T}{2N}\frac{T}{N}\sqrt{2}N^{3/2}=\mathcal{O}(\sqrt{\Delta t})$$ Which is the exponent you are obtaining. $$p=0.5$$ Keep in mind that the pointwise error, which must be of the order of the time step, is not the same as what you have computed, which is another error.

Answer 2

I cannot tell you whether that is the cause for your diminished order, but it is asking for trouble to add other approximations than just the one you want to do on purpose. In the current context, since you know $C_g$, there is no need to replace $C_g'$ by a finite difference quotient -- just use the explicit derivative.

Even better, introduce a new variable $\epsilon = EC_g$, express $\delta=\delta(E)=\delta(\epsilon/C_g)$ in terms of $\epsilon$, and then solve the ODE $$ \epsilon' = -\delta(\epsilon/C_g) $$ instead -- there is no ambiguity how it should be discretized.