Re: User supplied jacobian makes ode23s very slow

[Carlo De Falco]

On 6 Nov 2015, at 13:01, Alexander Bessman <address@hidden> wrote:

> Hello,
> 
> I'm using ode23s from odepkg to solve a simple 1D PDE, ∂u/∂t = ∂²u/dx².
> I have discretized it via the forward-time central-space scheme,
> yielding a tridiagonal system of ODEs, which ode23s has no problem
> solving. However, when I supply the system's Jacobian via odeset, the
> solution takes much longer to complete. The solution is still correct,
> but I would have expected it to be faster, not slower. I've included my
> code below. The result of running it is thus:
> 
> |Timesteps             |CPU-time            |Max stepsize        |
> |N  |with Jac |w/o Jac |with Jac  |w/o Jac  |with Jac  |w/o Jac  |
> |10 |432      |79      |1.3801 s  |1.1953 s |1.9914e-2 |2.0000e-1|
> |20 |1633     |98      |5.4225 s  |2.6037 s |5.2711e-3 |1.9447e-1|
> |40 |6442     |119     |26.7256 s |5.9159 s |8.4338e-4 |1.7272e-1|
> 
> I'd be grateful of someone could help me understand this phenomenon.
> //Alexander Bessman
> 
> -----------------------------------------------------------------------
> 
> pkg load odepkg
> 
> cpu = zeros(2,3);
> steps = zeros(2,3);
> htmax = zeros(2,3);
> iter = 0;
> for N = [10, 20, 40]    
>    % Jacobian
>    e = ones (N, 1);
>    jac = spdiags ([e, -2*e, e], [-1:1], N, N);
>    jac(end-1, end) = 2;
>    conf = odeset ('Jacobian', jac);
> 
>    u0 = zeros(N,1);
>    iter = iter + 1;
> 
>    tic
>    [t, y] = ode23s (@heat_eq, [0, 2], u0);
>    cpu(1, iter) = toc;
>    steps(1, iter) = length(t);
>    htmax(1, iter) = max(diff(t));
> 
>    tic
>    [t, y] = ode23s (@heat_eq, [0, 2], u0, conf);
>    cpu(2, iter) = toc;
>    steps(2, iter) = length(t);
>    htmax(2, iter) = max(diff(t));
> end
> 
> function dudt = heat_eq (t, u, kappa=1)
>    N = length(u);
>    e = ones (N, 1);
>    A = spdiags ([e, -2*e, e], [-1:1], N, N);
> 
>    % Boundary conditions
>    A(end-1, end) = 2;
>    b = zeros(N, 1);
>    if t <= 1
>        b(1) = 1;
>    end
> 
>    dudt = kappa * N^2 * (A*u + b);
> 
> 
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://lists.gnu.org/mailman/listinfo/help-octave



You have an incorrect jacobian which causes an unnecessary high number of 
iterations.
You have to multiply jac by N^2 * kappa, by doing so I get:

steps =

    79    98   119
    79    98   119

cpu =

   1.02229   2.04385   4.87316
   0.20588   0.24211   0.29443

cpu (1,:) ./cpu (2,:)
ans =

    4.9654    8.4419   16.5513

so solution with jacobian is up to 16x faster.
For larger problems the gain is even more impressive.

c.

ode23test.m
Description: ode23test.m

ATT00001.txt
Description: ATT00001.txt