[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: User supplied jacobian makes ode23s very slow
From: |
Carlo De Falco |
Subject: |
Re: User supplied jacobian makes ode23s very slow |
Date: |
Fri, 6 Nov 2015 16:48:21 +0000 |
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