Creating PDEs from ODEs

Contents

Creating PDEs from ODEs#

Modeling the Brusselator PDE:

\[\begin{split} \begin{align} \frac{\partial u}{\partial t} &= 1 + u^2v - 4.4u + \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) + f(x, y, t) \\ \frac{\partial v}{\partial t} &= 3.4u - u^2 v + \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) \end{align} \end{split}\]

import ComponentArrays.ComponentArray as CA
using SimpleUnPack
using OrdinaryDiffEq
using Plots

Model parameters and initial conditions

N=26
x_min=0
y_min=0
t_min=0
x_max=1
y_max=1
t_max=11.5
α=10.0
u0 = CA(u=zeros(N, N), v=zeros(N, N))
xx = range(x_min, x_max, length=N)
yy = range(y_min, y_max, length=N)
dx = step(xx)
dy = step(yy)
for i in 1:N, j in 1:N
    x = xx[j]
    y = yy[i]
    u0.u[i, j] = 22 * (y * (1 - y))
    u0.v[i, j] = 27 * (x * (1 - x))
end
u0

ComponentVector{Float64}(u = [0.0 0.0 … 0.0 0.0; 0.8447999999999999 0.8447999999999999 … 0.8447999999999999 0.8447999999999999; … ; 0.8448000000000007 0.8448000000000007 … 0.8448000000000007 0.8448000000000007; 0.0 0.0 … 0.0 0.0], v = [0.0 1.0368 … 1.0368000000000008 0.0; 0.0 1.0368 … 1.0368000000000008 0.0; … ; 0.0 1.0368 … 1.0368000000000008 0.0; 0.0 1.0368 … 1.0368000000000008 0.0])

Discretized PDE problem is a coupled ODE problem under the hood.

function model!(ds, s, p, t)
    @unpack α, xx, yy, dx, dy = p
    @unpack u, v = s
    brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5

    for (i, y) in enumerate(yy), (j, x) in enumerate(xx)
        uterm = 1.0 + v[i, j] * u[i, j]^2 - 4.4 * u[i, j] + brusselator_f(x, y, t)
        vterm = 3.4 * u[i, j] - v[i, j] * u[i, j]^2
        # No flux boundary conditions
        i_prev = clamp(i - 1, 1, N)
        i_next = clamp(i + 1, 1, N)
        j_prev = clamp(j - 1, 1, N)
        j_next = clamp(j + 1, 1, N)
        unabla = (u[i, j_prev] - 2u[i, j] + u[i, j_next]) / dx^2 + (u[i_prev, j] - 2u[i, j] + u[i_next, j]) / dy^2
        vnabla = (v[i, j_prev] - 2v[i, j] + v[i, j_next]) / dx^2 + (v[i_prev, j] - 2v[i, j] + v[i_next, j]) / dy^2
        ds.u[i, j] = α * unabla + uterm
        ds.v[i, j] = α * vnabla + vterm
    end
    return nothing
end

model! (generic function with 1 method)

ps = (;α, xx, yy, dx, dy)
tspan = (t_min, t_max)
prob = ODEProblem(model!, u0, tspan, ps)
@time sol = solve(prob, FBDF(), saveat=0.1)

 17.871480 seconds (56.04 M allocations: 2.687 GiB, 2.74% gc time, 44.19% compilation time)

retcode: Success
Interpolation: 1st order linear
t: 116-element Vector{Float64}:
  0.0
  0.1
  0.2
  0.3
  0.4
  0.5
  0.6
  0.7
  0.8
  0.9
  ⋮
 10.7
 10.8
 10.9
 11.0
 11.1
 11.2
 11.3
 11.4
 11.5
u: 116-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(u = ViewAxis(1:676, ShapedAxis((26, 26))), v = ViewAxis(677:1352, ShapedAxis((26, 26))))}}}}:
 ComponentVector{Float64}(u = [0.0 0.0 … 0.0 0.0; 0.8447999999999999 0.8447999999999999 … 0.8447999999999999 0.8447999999999999; … ; 0.8448000000000007 0.8448000000000007 … 0.8448000000000007 0.8448000000000007; 0.0 0.0 … 0.0 0.0], v = [0.0 1.0368 … 1.0368000000000008 0.0; 0.0 1.0368 … 1.0368000000000008 0.0; … ; 0.0 1.0368 … 1.0368000000000008 0.0; 0.0 1.0368 … 1.0368000000000008 0.0])
 ComponentVector{Float64}(u = [6.689207073770846 6.689207067610551 … 6.689207067603053 6.689207073763251; 6.6892070724373065 6.689207065259384 … 6.689207065251891 6.68920707242972; … ; 6.689207072432301 6.689207065254326 … 6.689207065246801 6.689207072424728; 6.689207073765753 6.689207067605476 … 6.689207067597996 6.689207073758142], v = [0.6994326397929781 0.699432646397615 … 0.6994326463996392 0.6994326397950068; 0.6994326414107599 0.6994326490329847 … 0.6994326490349373 0.6994326414128146; … ; 0.69943264141315 0.699432649035425 … 0.6994326490374251 0.6994326414151878; 0.6994326397954209 0.6994326464000422 … 0.6994326464020485 0.699432639797457])
 ComponentVector{Float64}(u = [6.301972198251998 6.301972199336778 … 6.301972199336681 6.301972198252081; 6.3019721980866095 6.301972199172716 … 6.301972199172935 6.301972198086271; … ; 6.301972198086817 6.301972199172501 … 6.301972199172979 6.301972198086416; 6.301972198251944 6.30197219933697 … 6.301972199336777 6.301972198252093], v = [0.5315881496698343 0.5315881482619097 … 0.5315881482619453 0.5315881496699136; 0.5315881495904892 0.5315881481807861 … 0.5315881481809179 0.5315881495905377; … ; 0.5315881495903905 0.5315881481809459 … 0.531588148180981 0.5315881495904846; 0.5315881496698469 0.531588148261821 … 0.5315881482619006 0.5315881496699069])
 ComponentVector{Float64}(u = [5.755196674136048 5.7551966740679745 … 5.755196674065452 5.755196674134518; 5.755196674152787 5.755196674083202 … 5.755196674085307 5.755196674153499; … ; 5.755196674154123 5.7551966740820895 … 5.755196674083953 5.7551966741553375; 5.755196674135575 5.755196674066988 … 5.755196674068343 5.755196674141106], v = [0.5757942156245692 0.575794215714415 … 0.5757942157147058 0.5757942156248221; 0.5757942156232319 0.5757942157139581 … 0.5757942157135367 0.5757942156234049; … ; 0.5757942156237544 0.5757942157139484 … 0.5757942157139765 0.5757942156230758; 0.5757942156255911 0.5757942157147256 … 0.5757942157142711 0.5757942156244138])
 ComponentVector{Float64}(u = [5.253187621170368 5.253187621177684 … 5.253187621183829 5.253187621176028; 5.253187621172077 5.2531876211773545 … 5.253187621180475 5.253187621173035; … ; 5.253187621175508 5.253187621180753 … 5.253187621177865 5.25318762116757; 5.253187621179091 5.253187621185167 … 5.2531876211775534 5.253187621164426], v = [0.627686075533669 0.627686075522572 … 0.6276860755219755 0.6276860755335099; 0.6276860755321874 0.6276860755207352 … 0.6276860755223885 0.6276860755318724; … ; 0.6276860755310556 0.6276860755211779 … 0.6276860755225785 0.6276860755330969; 0.6276860755323508 0.6276860755211132 … 0.6276860755234835 0.6276860755346332])
 ComponentVector{Float64}(u = [4.7948839562919945 4.794883956272841 … 4.794883956268465 4.7948839562888566; 4.794883956294568 4.794883956274686 … 4.794883956274557 4.794883956291749; … ; 4.794883956292303 4.794883956274119 … 4.794883956274626 4.794883956291318; 4.794883956288424 4.794883956271401 … 4.7948839562668315 4.794883956287106], v = [0.6838297779257222 0.6838297779494772 … 0.6838297779516055 0.6838297779273419; 0.683829777926216 0.6838297779494585 … 0.6838297779509365 0.6838297779279006; … ; 0.6838297779266246 0.683829777950256 … 0.6838297779520046 0.683829777928276; 0.6838297779259582 0.6838297779508118 … 0.6838297779518189 0.6838297779280206])
 ComponentVector{Float64}(u = [4.37624566220548 4.376245662203105 … 4.376245662187876 4.376245662195261; 4.376245662205875 4.37624566220286 … 4.376245662190342 4.376245662194081; … ; 4.376245662203153 4.376245662199941 … 4.376245662188757 4.376245662188912; 4.376245662201784 4.376245662199173 … 4.376245662185756 4.376245662188553], v = [0.7442956114202923 0.7442956114256083 … 0.7442956114441827 0.7442956114407644; 0.7442956114205009 0.744295611425057 … 0.7442956114422447 0.7442956114397004; … ; 0.7442956114305366 0.7442956114340723 … 0.7442956114430912 0.7442956114380267; 0.7442956114302246 0.7442956114336376 … 0.7442956114425566 0.7442956114384391])
 ComponentVector{Float64}(u = [3.9932930414747556 3.9932930414666976 … 3.99329304147827 3.993293041477775; 3.9932930414745704 3.9932930414663765 … 3.9932930414755226 3.9932930414761487; … ; 3.9932930414771772 3.993293041468205 … 3.99329304148653 3.9932930414902668; 3.993293041481713 3.99329304147128 … 3.9932930414896264 3.993293041489], v = [0.809240865534407 0.8092408655387715 … 0.8092408655106488 0.8092408655033402; 0.8092408655383015 0.8092408655431251 … 0.8092408655168128 0.8092408655075501; … ; 0.8092408655163589 0.809240865526568 … 0.8092408655118679 0.8092408655071202; 0.8092408655153632 0.8092408655220595 … 0.8092408655093597 0.809240865504734])
 ComponentVector{Float64}(u = [3.6421724511930913 3.6421724511907296 … 3.6421724511994378 3.6421724511960814; 3.6421724511892206 3.642172451196446 … 3.642172451198501 3.642172451196463; … ; 3.642172451183243 3.642172451193062 … 3.6421724511939684 3.6421724511967395; 3.642172451186445 3.64217245119141 … 3.642172451196388 3.642172451199673], v = [0.8790178643304916 0.8790178643289828 … 0.8790178643280716 0.8790178643313282; 0.8790178643311838 0.879017864324573 … 0.8790178643264001 0.8790178643326525; … ; 0.879017864340877 0.8790178643322286 … 0.8790178643167585 0.8790178643210865; 0.8790178643381278 0.8790178643366666 … 0.8790178643200341 0.8790178643203453])
 ComponentVector{Float64}(u = [3.3199291707012444 3.3199291707032548 … 3.3199291707010663 3.3199291707004113; 3.319929170701972 3.319929170702001 … 3.319929170700797 3.319929170701048; … ; 3.319929170697076 3.3199291706973217 … 3.3199291707045737 3.3199291707045284; 3.3199291706949436 3.319929170696432 … 3.3199291707043086 3.3199291707039325], v = [0.9539024882103184 0.9539024882090357 … 0.9539024882050198 0.9539024882056768; 0.9539024882093492 0.9539024882095543 … 0.9539024882047883 0.9539024882046976; … ; 0.9539024882113729 0.9539024882110108 … 0.9539024881986833 0.9539024881987345; 0.9539024882125682 0.953902488211188 … 0.9539024881985638 0.9539024881992971])
 ⋮
 ComponentVector{Float64}(u = [0.4051886081187576 0.4051791710730025 … 0.40334891175645804 0.4033342431831844; 0.405220476980117 0.405211246578527 … 0.4033565403218058 0.40334168466845643; … ; 0.4070182809267365 0.4070146890105795 … 0.4040155013419711 0.4039917199605697; 0.4069904513550748 0.4069853608305007 … 0.40401675854132785 0.40399312670480125], v = [4.117488331862 4.117488316969156 … 4.117486071495318 4.117486053790609; 4.117488347758788 4.11748833287931 … 4.117486079977186 4.117486062198347; … ; 4.117489713559648 4.117489697801436 … 4.117487019304264 4.117486997546294; 4.117489715656763 4.117489699792559 … 4.117487024520878 4.117487002794964])
 ComponentVector{Float64}(u = [0.4089050031558753 0.4088955485799748 … 0.4070626698456526 0.4070479810540002; 0.40893689427320534 0.40892764635357814 … 0.4070703081219203 0.4070554321417236; … ; 0.4107365503001098 0.4107329385555746 … 0.40773033031316913 0.4077065222290802; 0.4107087206313452 0.41070361006877204 … 0.4077315931688002 0.4077079346716665], v = [4.187086937592285 4.187086940462689 … 4.187087346872905 4.187087349630567; 4.187086931063636 4.187086933937803 … 4.187087345531992 4.187087348324864; … ; 4.187086426069602 4.187086430400608 … 4.18708720951285 4.1870872147391465; 4.1870864281643785 4.187086432598891 … 4.187087208976772 4.187087214192311])
 ComponentVector{Float64}(u = [0.41337159552675873 0.41336212202713746 … 0.41152641635972304 0.4115117057403307; 0.41340351065993214 0.41339424382389844 … 0.41153406509608587 0.4115191671783728; … ; 0.41520516488575965 0.41520153174740526 … 0.4121952326969398 0.4121713958063902; 0.4151773350959456 0.4151722029103023 … 0.41219650167096133 0.41217281440499487], v = [4.2555166605243215 4.255516682569915 … 4.2555199508120545 4.255519975659536; 4.255516629797725 4.255516651840962 … 4.2555199388906635 4.255519963884101; … ; 4.255514106263738 4.255514132270548 … 4.255518642311936 4.255518676649856; 4.255514108369646 4.255514134710031 … 4.25551863555573 4.255518669843348])
 ComponentVector{Float64}(u = [0.41858558805485857 0.4185760942482331 … 0.4167373589901289 0.41672262498726265; 0.41861752889796067 0.4186082417518908 … 0.41674501893994126 0.4167300975216002; … ; 0.42042132268356447 0.42041766657883056 … 0.4174074144510824 0.41738354671889616; 0.42039349283755917 0.42038833743430193 … 0.41740869001125847 0.4173849719454183], v = [4.322721160462818 4.3227212030859965 … 4.322727538471145 4.3227275869855815; 4.322721103832377 4.322721146454261 … 4.322727515203276 4.322727563980549; … ; 4.322716418988498 4.322716468263698 … 4.3227249743898355 4.322725039893186; 4.322716421033813 4.322716470893758 … 4.322724960938186 4.3227250263524315])
 ComponentVector{Float64}(u = [0.4245329670296315 0.4245234514828522 … 0.4226814740780798 0.4226667150615438; 0.42456493535134865 0.42455562647294587 … 0.42268914603280217 0.4226741994757777; … ; 0.4263710176320815 0.4263673369527163 … 0.4233528558374532 0.4233289551172411; 0.4263431877401922 0.4263380075051826 … 0.42335413844965447 0.42333038742659845], v = [4.388630465243695 4.388630529895587 … 4.388640147598591 4.388640221429119; 4.388630380926063 4.388630445569856 … 4.388640112184344 4.3886401864043565; … ; 4.3886233843368885 4.388623458510161 … 4.388636239596253 4.388636338433819; 4.388623386303063 4.388623461317936 … 4.388636218973526 4.388636317694105])
 ComponentVector{Float64}(u = [0.4312348611394853 0.431225322357606 … 0.42937987584921516 0.42936509005063084; 0.43126685883001037 0.43125752672501333 … 0.42938756063683264 0.42937258716310517; … ; 0.4330753888968429 0.4330716819540236 … 0.4300526769519466 0.43002874092560356; 0.4330475589557585 0.4330423521774892 … 0.43005396709339583 0.4300301808116421], v = [4.453148213747968 4.453148301943742 … 4.453161431606219 4.453161532540756; 4.453148099838494 4.453148188018988 … 4.4531613832110395 4.453161484670716; … ; 4.453138630549306 4.453138731331823 … 4.453156085488993 4.453156220006648; 4.453138632430058 4.453138734334866 … 4.45315605721192 4.453156191563096])
 ComponentVector{Float64}(u = [0.43878231712748117 0.43877275347247613 … 0.4369235916210962 0.4369087771246495; 0.438814346215227 0.4388049892273468 … 0.4369312901452966 0.43691628782553177; … ; 0.44062549547409874 0.4406217603875273 … 0.43759791329094316 0.43757393945852197; 0.4405976655761481 0.4405924303348423 … 0.4375992115007114 0.43757538749849895], v = [4.516071689112847 4.516071802515063 … 4.516088693556699 4.516088823536312; 4.516071543569633 4.516071656962706 … 4.516088631265929 4.51608876191583; … ; 4.5160594285280675 4.516059557823254 … 4.516081806647057 4.516081979371033; 4.516059430223459 4.516059560961294 … 4.516081770169277 4.516081942638717])
 ComponentVector{Float64}(u = [0.4472637456926849 0.44725415521708206 … 0.4454009936251712 0.44538614826662054; 0.447295808552452 0.447286424752752 … 0.4454087069466673 0.4453936736105617; … ; 0.44910977614176506 0.4491060107399681 … 0.4460769525409123 0.44605293805454355; 0.4490819462919651 0.4490766804205404 … 0.446078259465981 0.44605439485538284], v = [4.577152756980482 4.577152897559382 … 4.577173837831207 4.577173999044842; 4.5771525774087864 4.577152717971668 … 4.577173760571403 4.5771739226112675; … ; 4.577137615424292 4.577137775432915 … 4.577165291943445 4.577165505749394; 4.577137616917139 4.577137778683609 … 4.577165246603297 4.5771654601096285])
 ComponentVector{Float64}(u = [0.4566119673191592 0.4566023482443103 … 0.4547449214035227 0.4547300431301215; 0.4566440661887783 0.45663465380312057 … 0.45475265050501895 0.45473758409092446; … ; 0.458461038571534 0.45845724086177464 … 0.45542262643751885 0.45539856860445926; 0.4584332087603456 0.45842791024067503 … 0.4554239426546255 0.4554000347426254], v = [4.6365793658403724 4.636579535396861 … 4.636604793529522 4.636604988053902; 4.636579149984905 4.636579319513288 … 4.636604700304777 4.636604895820478; … ; 4.63656115258759 4.6365613453292385 … 4.636594478511001 4.636594736120983; 4.636561153876165 4.636561348716321 … 4.636594423725848 4.636594680991979])

Runtime environment#

using InteractiveUtils
InteractiveUtils.versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 4 default, 1 interactive, 4 GC (on 4 virtual cores)
Environment:
  JULIA_CPU_TARGET = generic;icelake-server,clone_all;znver3,clone_all
  JULIA_CONDAPKG_OFFLINE = true
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_CI = true
  LD_LIBRARY_PATH = /opt/hostedtoolcache/Python/3.13.9/x64/lib
  JULIA_NUM_THREADS = auto

using Pkg
Pkg.status()

Status `~/work/jl-pde/jl-pde/Project.toml`
  [b0b7db55] ComponentArrays v0.15.30
⌃ [aae7a2af] DiffEqFlux v4.4.0
  [5b8099bc] DomainSets v0.7.16
  [7da242da] Enzyme v0.13.108
  [f6369f11] ForwardDiff v1.3.0
  [d3d80556] LineSearches v7.5.1
⌃ [b2108857] Lux v1.21.1
  [94925ecb] MethodOfLines v0.11.9
  [961ee093] ModelingToolkit v10.31.0
  [315f7962] NeuralPDE v5.20.0
  [8913a72c] NonlinearSolve v4.12.0
⌅ [7f7a1694] Optimization v4.8.0
⌃ [36348300] OptimizationOptimJL v0.4.5
⌃ [42dfb2eb] OptimizationOptimisers v0.3.11
⌃ [500b13db] OptimizationPolyalgorithms v0.3.1
  [1dea7af3] OrdinaryDiffEq v6.103.0
  [91a5bcdd] Plots v1.41.2
  [ce78b400] SimpleUnPack v1.1.0
  [37e2e46d] LinearAlgebra v1.12.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This notebook was generated using Literate.jl.