Solving Partial Differential Equations with Julia
August 8 2018 in Differential Equations, HPC, Julia, Mathematics, Programming, Uncategorized | Tags: | Author: Christopher Rackauckas
Here is a talk from JuliaCon 2018 where I describe how to use the tooling across the Julia ecosystem to solve partial differential equations (PDEs), and how the different areas of the ecosystem are evolving to give top-notch PDE solver support.
Sandeep Reddy Bukka
says:Hi Chris
Thanks for the reply.
I was able to figure out the error. I think it is related to the new syntax of julia. I think i missed a broadcast equal to.
I was using du = L*u, instead of du.=L*u.
It works now !
Regards
Sandeep
Sandeep Reddy Bukka
says:Hi Christopher
Great work
The video is amazing. However i would like to mention that when i try to replicate the example of spectral time stepping of the heat equation using the current version of julia 1.4.2. I am getting some errors. I have changed the old commands of A_mul_B! but still i am getting errors in the solver convergence. Specifically, i am getting this
[CVODES ERROR] CVode
At t = 0 and h = 1.3093e-009, the corrector convergence test failed repeatedly or with |h| = hmin.
Can you please update on this?. It would be of great help.
Thanks
Sandeep
Christopher Rackauckas
says:My guess is that something else changed as well. The video is quite old. Open an issue or post some code into the Julia Discourse and we can work this out.
Sean
says:Is there perhaps a tutorial or example of this somewhere (especially for the case of non-linear reaction terms)?
Christopher Rackauckas
says:Not quite yet. I want to write one this weekend though.
Sean
says:Hi Chris
Nice talk!
Connecting with what you said at the end of your talk, I have seen that there are now wrappers for 2D and 3D FEM heat equation problems (http://docs.juliadiffeq.org/v1.10.1/tutorials/femheat_example.html) but is there a wrapper for the (1D) multi-variable (i.e. reaction-diffusion) case?
Christopher Rackauckas
says:Hey, that’s actually a very old version of the documentation (from way before this talk!). We’ve gone back to the drawing board and are taking a more modular approach. Instead, what we recommend is discretizing using tools like DiffEqOperators.jl (https://github.com/JuliaDiffEq/DiffEqOperators.jl) or ApproxFun.jl (https://github.com/JuliaApproximation/ApproxFun.jl), and then solving the resulting ODEs/nonlinear equations/linear equations.
Hesham M
says:Interesting video. I’m still going through it.
In block [46] in the jupyter notebook, I think you had the order and derivative swapped:
A = DerivativeOperator{Float64}(order,deriv,Δx,N,:Dirichlet0,:Dirichlet0)
Thanks for sharing!
Yoshitaka Taguchi
says:Hi chris, I watched your talk on youtube. Thank you for great lecture!
In the talk, you showed that some PDEs come down to ODE, like u_t = u_xx + f(u, t). Converting second derivetive of x to a matrix, this PDE becomes ODE, whose variable is t.
However, we can take x as variable of ODE. For example, let v = u_x. Then this PDE become simultaneous ODE of u_x = v and v_x = u_t – f(u,t).
If the range of x is very big and range of t is small, I thought taking x as variable of ODE is better because we don’t need to prepare big matrix for differentiation of x.
What do you think about this method? Is this approach efficient?
Christopher Rackauckas
says:Your suggestion ends up solving a slightly different problem because it requires a different initial condition. For your case, your initial condition would need to be u(0,t), that is the value at x=0 for all time points. The discretization I showed is for when you have u(x,0), the value at t=0. Normally this is the initial condition that the problem is given with.
Dimiter Georgiev
says:Hi Chriss, great video!
from your API i see first and second order operarator can be construccted as of eg.
D2 = DerivativeOperator{Float64}(2,2,1/100,11,:Neumann,:Dirichlet; BC=(0.0,0.0))
D1 = DerivativeOperator{Float64}(1,2,1/100,11,:Neumann,:Dirichlet; BC=(0.0,0.0))
but how should i construct something like
L = D2 + 1/x * D1
with boundaries?
I expect somehow first L has to be defined and then some function
to construct the boundaries be applied?
Christopher Rackauckas
says:Note that in the video I mention that it’s not entirely ready for release. That’s because we haven’t completely worked out the composition stuff yet.
fkoessel
says:I’ve watched your talk on YouTube. It was a comprehensible, well structured, brilliant talk. I enjoyed it very much!
Are your slides online by any chance?
Good work!
Christopher Rackauckas
says:Here you go: http://nbviewer.jupyter.org/github/ChrisRackauckas/JuliaCon2018PDEWorkshop/blob/master/PDEWorkshop.ipynb