search

SciML / PDERoadmap

A repository for the discussion of PDE tooling for scientific machine learning (SciML) and physics-informed machine learning
https://tutorials.sciml.ai/
18 stars 6 forks source link

Test variations on DifferentialEquations parameters for the PDE proof of concept #25

Closed jlperla closed 6 years ago

jlperla commented 6 years ago

The "stability" question with the dynamics is whether we can go backwards in time for significant $t$ and not have weird things cropping up. Plotting with some of the solvers you can see that you start losing the ordinal order, etc.

When the base code in #23 is done, we should try the following to see if stability increases, etc:

jlperla commented 6 years ago

The AutoTsit5 and Rosenbrock23() are unstable for large time-domains (whatever the reason) which you can see with the current settings in the transition_dynamics_example.jl Not sure why the ImplicitEuler works better, but it seems to for the "not-time-varying" case at least.

jlperla commented 6 years ago

Let me take that back. Rosenbrock23() seems to be stable in the non-time varying case when the adaptive grid is taken off (with a dt =0.01). It could just be that the adaptive time-grid is the issue rather than the algorithm itself?

ChrisRackauckas commented 6 years ago

The test case isn't a good test case for a normal integration. You're testing how what a solver does at a steady state. What's happening is that it sees that the derivative is zero, so it just exponentially grows the stepsize until numerical error dominates and it falls off. ImplicitEuler() is an interesting case because when the implicit part converges it will essentially ensure that it has to not be changing, so it's not going to accumulate numerical error. But any other algorithm (even other SDIRKs) are going to be subject to this problem unless there's a dtmax which is set. But this is just some bizarre behavior of adaptivity if you sit at steady states.

CVODE_BDF() doesn't display this behavior, so I wonder if they are doing something special to cancel out that numerical error or if it's just a property of BDF. We'll see when @YingboMa 's GSoC is done. Probably just needs a tweak to something small like the derivative calculation.

For backwards integration just reverse tspan. For example:

using DifferentialEquations

function time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
    M_bar = M + 2
    x = linspace(x_min, x_max, M)
    Delta = x[2]-x[1]
    mu_tilde_T = -0.1+T+0.1*x
    sigma_tilde_T = sigma_bar*x
    c_tilde = exp.(x)

    L_1_plus = Tridiagonal(zeros(M_bar-1), -1 * ones(M_bar), ones(M_bar-1))[2: end-1, :]/Delta
    L_1_plus[:, 2] = L_1_plus[:, 2] + L_1_plus[:, 1]
    L_1_plus[:, end-1] = L_1_plus[:, end-1] + L_1_plus[:, end]
    L_1_plus = L_1_plus[:, 2: end-1]

    L_2 = Tridiagonal(ones(M_bar-1), -2 * ones(M_bar), ones(M_bar-1))[2: end-1, :]/(Delta^2)
    L_2[:, 2] = L_2[:, 2] + L_2[:, 1]
    L_2[:, end-1] = L_2[:, end-1] + L_2[:, end]
    L_2 = L_2[:, 2: end-1]

    L_tilde = Diagonal(rho*ones(M)) - (Diagonal(mu_tilde_T)*L_1_plus + Diagonal(sigma_tilde_T.^2/2)*L_2)

    u_T = L_tilde\c_tilde

    function f(du,u,p,t)
        mu_tilde = mu_tilde_T
        L = Diagonal(rho*ones(M)) - (Diagonal(mu_tilde)*L_1_plus + Diagonal(sigma_tilde_T.^2/2)*L_2)
        A_mul_B!(du,L,u)
        du .+= c_tilde
    end

    tstops=(T , 0.0)
    prob = ODEProblem(f, u_T, tstops)
    solve(prob, algorithm)
end

x_min = 0.01
x_max = 1
M = 40
T = 100.0
N = 10 #This is not being used right now.
rho = 0.05
sigma_bar = 0.1
algorithm = ImplicitEuler()

sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)

Notice that I also changed the PDE so that way it's doing + on the constant. When that's done, the reverse PDE has a steady state and it stays there.

jlperla commented 6 years ago

Perfect, thanks. You beat me to asking you the question on the divergence, and your explanation makes sense. Steven, can you tell me when you have your updates to the operators pushed to the server so I can play around with it?

ChrisRackauckas commented 6 years ago

FWIW the script I was testing the different methods with was:

using DifferentialEquations

function time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
    M_bar = M + 2
    x = linspace(x_min, x_max, M)
    Delta = x[2]-x[1]
    mu_tilde_T = -0.1+T+0.1*x
    sigma_tilde_T = sigma_bar*x
    c_tilde = exp.(x)

    L_1_plus = Tridiagonal(zeros(M_bar-1), -1 * ones(M_bar), ones(M_bar-1))[2: end-1, :]/Delta
    L_1_plus[:, 2] = L_1_plus[:, 2] + L_1_plus[:, 1]
    L_1_plus[:, end-1] = L_1_plus[:, end-1] + L_1_plus[:, end]
    L_1_plus = L_1_plus[:, 2: end-1]

    L_2 = Tridiagonal(ones(M_bar-1), -2 * ones(M_bar), ones(M_bar-1))[2: end-1, :]/(Delta^2)
    L_2[:, 2] = L_2[:, 2] + L_2[:, 1]
    L_2[:, end-1] = L_2[:, end-1] + L_2[:, end]
    L_2 = L_2[:, 2: end-1]

    L_tilde = Diagonal(rho*ones(M)) - (Diagonal(mu_tilde_T)*L_1_plus + Diagonal(sigma_tilde_T.^2/2)*L_2)

    u_T = L_tilde\c_tilde

    function f(du,u,p,t)
        #mu_tilde = -0.1+0.01*t+0.1*x#mu_tilde_T#-0.1+t+0.1*x
        mu_tilde = mu_tilde_T#-0.1+t+0.1*x
        L = Diagonal(rho*ones(M)) - (Diagonal(mu_tilde)*L_1_plus + Diagonal(sigma_tilde_T.^2/2)*L_2)
        A_mul_B!(du,L,u)
        du .-= c_tilde
    end

    tstops=(0.0, T)
    prob = ODEProblem(f, u_T, tstops)
    solve(prob, algorithm)
end

x_min = 0.01
x_max = 1
M = 40
T = 100.0
N = 10 #This is not being used right now.
rho = 0.05
sigma_bar = 0.1
#algorithm = Tsit5()
#algorithm = AutoTsit5(Rosenbrock23())

algorithm = ImplicitEuler()
sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
plot(sol, vars=1:5:M)

algorithm = Rosenbrock23()
sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
plot(sol, vars=1:5:M)

algorithm = Rodas4()
sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
plot(sol, vars=1:5:M)

algorithm = KenCarp4()
sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
plot(sol, vars=1:5:M)

algorithm = TRBDF2()
sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
plot(sol, vars=1:5:M)

algorithm = CVODE_BDF()
sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar, algorithm)
plot(sol, vars=1:5:M)

The fact that it happens both on Rosenbrock and SDIRK methods tells me it's not something with implicit stepping handling since Rosenbrock isn't implicit (just semi-implicit). It also shows that it's method dependent and fixable since Sundials doesn't display this behavior. It's highly related to some behaviors seen on this benchmark problem:

http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/ROBER.ipynb

So yeah, it's now you it's me. It's on our end to work out a more numerically stable way of handling this. @YingboMa should be aware of this though I am not sure of a solution quite yet. But as mentioned earlier it's a behavior that happens as the derivative term approaches zero and thus the error -> 0 causing time steps to go infinitely large, so it would be nice to fix but it's constrained to somewhat of a corner case.

Another interesting fact is that both finite and autodiff fall off at the same point. My guess is that it's a side effect of PI timestepping adaptivity.

jlperla commented 6 years ago

@ChrisRackauckas Thanks, sounds great. I don't think this is a big issue in the short term because we will only really be solving problems with sufficient time variation at this point. But worth noting for diagonstics, if nothing else.

jlperla commented 6 years ago

One comment on this "bug". The stability issues (even for a very large T) may not be cropping up now that we have stopped swapping the sign of the operator. The latest in https://github.com/JuliaDiffEq/PDERoadmap/blob/master/transition_dynamics_examples/transition_dynamics_example.jl seems very stable for all sorts of parameter and time variation (or lack thereof) and the signs of everything make sense.