[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: User supplied jacobian makes ode23s very slow
From: |
Alexander Bessman |
Subject: |
Re: User supplied jacobian makes ode23s very slow |
Date: |
Sat, 07 Nov 2015 15:13:06 +0100 |
On Fri, 2015-11-06 at 16:48 +0000, Carlo De Falco wrote:
>
> 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.
>
Thank you. Such a silly mistake, but I suppose I should have expected as
much. Very impressive speed gains indeed!
//Alexander