From f0885a87add238185f6db50ba48608751546dc07 Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Sun, 3 Aug 2025 19:50:59 +0200 Subject: [PATCH 1/9] Add tests of IDEs (integro-diff-eq) with integrals [0..x] and [x..END] --- test/pde_systems/MOL_1D_Integration.jl | 131 +++++++++++++++++++++++++ 1 file changed, 131 insertions(+) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index aef8bab5d..65275aa08 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -159,3 +159,134 @@ end # @test cumuSumsol ≈ exact atol = 0.36 end + +@testset "Test 03: Partial integro-differential equation with time derivative, with sys transformation" begin + # Equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t)) + # Initial condition: u(0,x) = 0 + # Boundary condition: u(t,0) = 0 + # Analytical solution: u(t,x) = x * (1 - exp(-t)) + + @parameters t, x + @variables u(..) + Dt = Differential(t) + Dx = Differential(x) + xmin = 0.0 + xmax = 2.0 + + Ix_0_to_x = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x + + #eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t))] + eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t - 0.2))] + bcs = [ + # u(0,x) ~ 0.0, # Initial condition + u(0,x) ~ x * (1 - exp(-0.2)), # Initial condition + u(t,0) ~ 0.0] # Boundary condition + domains = [ + t ∈ Interval(0.0, 1.0), + x ∈ Interval(xmin, xmax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, x], [u(t, x)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([x => 50], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + xdisc = sol[x] + tdisc = sol[t] + usol = sol[u(t,x)] + + # Calculate analytical solution at grid points + analytical_solution(t, x) = x * (1 - exp(-(t+0.2)) ) + exact = [analytical_solution(t_, x_) for t_ in tdisc, x_ in xdisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.1 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(xdisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="x", + # ylabel="u(t,x)", + # title="Solution at t = $(round(t_val, digits=3))", + # ylim=(0, maximum(usol)*1.1), + # xlim=(xmin, xmax)) + + # plot!(xdisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + + # #legend!(:topright) + # end + # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation +end + +@testset "Test 03b: Mirrored version of: Partial integro-differential equation with time derivative, with sys transformation" begin + # Original equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t)) + # Mirrored version with transformation y = -x, x = -y: + # Equation: ∂u/∂t + u(t,y) + ∫₀ˣ u(t,ξ) dξ = -x + x²/2 * (1 - exp(-t)) + # Initial condition: u(0,y) = 0 + # Boundary condition: u(t,0) = 0 + # Analytical solution: u(t,y) = -y * (1 - exp(-t)) + + @parameters t, y + @variables u(..) + Dt = Differential(t) + Dy = Differential(y) + ymin = -2.0 + ymax = 0.0 + + # Integral from y to 0 (note the reversed limits) + Iy_to_0 = Integral(y in DomainSets.ClosedInterval(y, ymax)) + + #eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t))] + eqs = [Dt(u(t,y)) + u(t,y) - Iy_to_0(u(t,y)) ~ -y + y^2/2*(1-exp(-t))] + bcs = [ + u(0,y) ~ 0.0, # Initial condition + u(t,0) ~ 0.0] # Boundary condition (at y = 0, which corresponds to x = 0) + domains = [ + t ∈ Interval(0.0, 1.0), + y ∈ Interval(ymin, ymax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, y], [u(t, y)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([y => 50], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + ydisc = sol[y] + tdisc = sol[t] + usol = sol[u(t,y)] + + # Calculate analytical solution at grid points + # Since u(t,x) = x * (1 - exp(-t)) and y = -x, then u(t,y) = -y * (1 - exp(-t)) + analytical_solution(t, y) = -y * (1 - exp(-t)) + exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.1 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(ydisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="x", + # ylabel="u(t,x)", + # title="Solution at t = $(round(t_val, digits=3))", + # ylim=(0, maximum(usol)*1.1), + # xlim=(xmin, xmax)) + + # plot!(xdisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + # end + # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation +end From eb5d7e68216392bc10962f6cd21bb32d093c83ee Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Sun, 3 Aug 2025 22:22:25 +0200 Subject: [PATCH 2/9] Add tests of PIDEs (integro-pde) with integrals [0..x] and [x..END] --- test/pde_systems/MOL_1D_Integration.jl | 149 ++++++++++++++++++++++++- 1 file changed, 147 insertions(+), 2 deletions(-) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index 65275aa08..4f8d8b3c5 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -219,8 +219,7 @@ end # label="Analytical", # linewidth=2, # linestyle=:dash) - - # #legend!(:topright) + # end # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation end @@ -290,3 +289,149 @@ end # end # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation end + +@testset "Test 04: PIDE with spatial derivative and source term - domain [0, 4π]" begin + # Equation: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + # where h(x,t) = -sin(x)sin(t) + sin(x)cos(t) - cos(x)cos(t) + cos(t)(1 - cos(x)) + # On domain: x ∈ [0,L]=[0,2π] and t ∈ [0,T]=[0,1] + # Initial condition: u(0,x) = sin(x) + # Boundary conditions: u(t,0) = 0, u(t,2π) = sin(2π)cos(t) = 0 + # Analytical solution: u(t,x) = sin(x)cos(t) + + @parameters t, x + @variables u(..) + Dt = Differential(t) + Dx = Differential(x) + xmin = 0.0 + xmax = 4π + + I0_to_x = Integral(x in DomainSets.ClosedInterval(xmin, x)) # integral from 0 to x + + # Source term h(x,t) that admits the analytical solution + h(x, t) = -sin(x)*sin(t) + sin(x)*cos(t) - cos(x)*cos(t) + cos(t)*(1 - cos(x)) + + # PIDE: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + eqs = [Dt(u(t,x)) ~ Dx(u(t,x)) - u(t,x) - I0_to_x(u(t,x)) + h(x, t)] + bcs = [ + u(0, x) ~ sin(x), # Initial condition + u(t, 0) ~ 0.0, # Boundary condition at x = 0 + u(t, xmax) ~ sin(xmax)*cos(t) # Boundary condition at x = 4π (sin(4π) = 0) + ] + domains = [ + t ∈ Interval(0.0, 1.0), + x ∈ Interval(xmin, xmax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, x], [u(t, x)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([x => 200], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + xdisc = sol[x] + tdisc = sol[t] + usol = sol[u(t,x)] + + # Calculate analytical solution at grid points + analytical_solution(t, x) = sin(x) * cos(t) + exact = [analytical_solution(t_, x_) for t_ in tdisc, x_ in xdisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.36 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(xdisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="x", + # ylabel="u(t,x)", + # title="PIDE Solution [0,4π] at t = $(round(t_val, digits=3))", + # ylim=(-1.1, 1.1), + # xlim=(xmin, xmax)) + + # plot!(xdisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + + # end + # gif(anim, "pide_solution_4pi.gif", fps=10) # Save the animation +end + + +@testset "Test 04b: PIDE with spatial derivative and source term - mirrored domain [-4π, 0]" begin + # Mirrored version with transformation y = -x, so x = -y + # Original equation: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + # Transformed equation: ∂u/∂t = -∂u/∂y - u(t,y) - ∫ₓ⁰ u(t,η) dη + h(-y,t) + # where h(-y,t) = -sin(-y)sin(t) + sin(-y)cos(t) - cos(-y)cos(t) + cos(t)(1 - cos(-y)) + # = sin(y)sin(t) - sin(y)cos(t) - cos(y)cos(t) + cos(t)(1 - cos(y)) + # On domain: y ∈ [-4π, 0] and t ∈ [0,1] + # Initial condition: u(0,y) = sin(-y) = -sin(y) + # Boundary conditions: u(t,0) = 0, u(t,-4π) = sin(-(-4π))cos(t) = sin(4π)cos(t) = 0 + # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) + + @parameters t, y + @variables u(..) + Dt = Differential(t) + Dy = Differential(y) + ymin = -4π + ymax = 0.0 + + # Integral from y to 0 (note the reversed limits) + Iy_to_0 = Integral(y in DomainSets.ClosedInterval(y, ymax)) + + # Source term h(-y,t) transformed for y coordinate + h_transformed(y, t) = sin(y)*sin(t) - sin(y)*cos(t) - cos(y)*cos(t) + cos(t)*(1 - cos(y)) + + # PIDE: ∂u/∂t = -∂u/∂y - u(t,y) - ∫ʸ⁰ u(t,η) dη + h(-y,t) + eqs = [Dt(u(t,y)) ~ -Dy(u(t,y)) - u(t,y) - Iy_to_0(u(t,y)) + h_transformed(y, t)] + bcs = [ + u(0, y) ~ -sin(y), # Initial condition: u(0,y) = sin(-y) = -sin(y) + u(t, 0) ~ 0.0, # Boundary condition at y = 0 + u(t, ymin) ~ -sin(ymin)*cos(t) # Boundary condition at y = -4π (sin(-4π) = 0) + ] + domains = [ + t ∈ Interval(0.0, 1.0), + y ∈ Interval(ymin, ymax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, y], [u(t, y)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([y => 200], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + ydisc = sol[y] + tdisc = sol[t] + usol = sol[u(t,y)] + + # Calculate analytical solution at grid points + # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) + analytical_solution(t, y) = -sin(y) * cos(t) + exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.36 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(ydisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="y", + # ylabel="u(t,y)", + # title="PIDE Solution [-4π,0] at t = $(round(t_val, digits=3))", + # ylim=(-1.1, 1.1), + # xlim=(ymin, ymax)) + + # plot!(ydisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + + # end + # gif(anim, "pide_solution_mirrored_4pi.gif", fps=10) # Save the animation +end From 19a0c7312529f91c6f4d1011cb4eaf0cf94573e3 Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Sun, 3 Aug 2025 22:32:38 +0200 Subject: [PATCH 3/9] Mark currently broken tests (Int[x..END]) as broken --- test/pde_systems/MOL_1D_Integration.jl | 130 +++++++++++++------------ 1 file changed, 66 insertions(+), 64 deletions(-) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index 4f8d8b3c5..1db4ecfdd 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -254,40 +254,41 @@ end # Discretize and compute numerical solution disc = MOLFiniteDifference([y => 50], t) - prob = discretize(pde_system, disc) - sol = solve(prob, Tsit5(), saveat=0.1) - - # Extract numerical solution at grid points - ydisc = sol[y] - tdisc = sol[t] - usol = sol[u(t,y)] - - # Calculate analytical solution at grid points - # Since u(t,x) = x * (1 - exp(-t)) and y = -x, then u(t,y) = -y * (1 - exp(-t)) - analytical_solution(t, y) = -y * (1 - exp(-t)) - exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] - - # Compare numerical and analytical solutions - @test usol ≈ exact atol=0.1 + @test_broken (discretize(pde_system, disc) isa ODEProblem) + # prob = discretize(pde_system, disc) + # sol = solve(prob, Tsit5(), saveat=0.1) - # # Create animated plot - # using Plots # Add this import - # anim = @animate for (i, t_val) in enumerate(tdisc) - # plot(ydisc, usol[i, :], - # label="Numerical", - # linewidth=2, - # xlabel="x", - # ylabel="u(t,x)", - # title="Solution at t = $(round(t_val, digits=3))", - # ylim=(0, maximum(usol)*1.1), - # xlim=(xmin, xmax)) + # # Extract numerical solution at grid points + # ydisc = sol[y] + # tdisc = sol[t] + # usol = sol[u(t,y)] + + # # Calculate analytical solution at grid points + # # Since u(t,x) = x * (1 - exp(-t)) and y = -x, then u(t,y) = -y * (1 - exp(-t)) + # analytical_solution(t, y) = -y * (1 - exp(-t)) + # exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # # Compare numerical and analytical solutions + # @test usol ≈ exact atol=0.1 + + # # # Create animated plot + # # using Plots # Add this import + # # anim = @animate for (i, t_val) in enumerate(tdisc) + # # plot(ydisc, usol[i, :], + # # label="Numerical", + # # linewidth=2, + # # xlabel="x", + # # ylabel="u(t,x)", + # # title="Solution at t = $(round(t_val, digits=3))", + # # ylim=(0, maximum(usol)*1.1), + # # xlim=(xmin, xmax)) - # plot!(xdisc, exact[i, :], - # label="Analytical", - # linewidth=2, - # linestyle=:dash) - # end - # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation + # # plot!(xdisc, exact[i, :], + # # label="Analytical", + # # linewidth=2, + # # linestyle=:dash) + # # end + # # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation end @testset "Test 04: PIDE with spatial derivative and source term - domain [0, 4π]" begin @@ -399,39 +400,40 @@ end # Discretize and compute numerical solution disc = MOLFiniteDifference([y => 200], t) - prob = discretize(pde_system, disc) - sol = solve(prob, Tsit5(), saveat=0.1) - - # Extract numerical solution at grid points - ydisc = sol[y] - tdisc = sol[t] - usol = sol[u(t,y)] - - # Calculate analytical solution at grid points - # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) - analytical_solution(t, y) = -sin(y) * cos(t) - exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] - - # Compare numerical and analytical solutions - @test usol ≈ exact atol=0.36 + @test_broken (discretize(pde_system, disc) isa ODEProblem) + # prob = discretize(pde_system, disc) + # sol = solve(prob, Tsit5(), saveat=0.1) - # # Create animated plot - # using Plots # Add this import - # anim = @animate for (i, t_val) in enumerate(tdisc) - # plot(ydisc, usol[i, :], - # label="Numerical", - # linewidth=2, - # xlabel="y", - # ylabel="u(t,y)", - # title="PIDE Solution [-4π,0] at t = $(round(t_val, digits=3))", - # ylim=(-1.1, 1.1), - # xlim=(ymin, ymax)) + # # Extract numerical solution at grid points + # ydisc = sol[y] + # tdisc = sol[t] + # usol = sol[u(t,y)] + + # # Calculate analytical solution at grid points + # # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) + # analytical_solution(t, y) = -sin(y) * cos(t) + # exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # # Compare numerical and analytical solutions + # @test usol ≈ exact atol=0.36 + + # # # Create animated plot + # # using Plots # Add this import + # # anim = @animate for (i, t_val) in enumerate(tdisc) + # # plot(ydisc, usol[i, :], + # # label="Numerical", + # # linewidth=2, + # # xlabel="y", + # # ylabel="u(t,y)", + # # title="PIDE Solution [-4π,0] at t = $(round(t_val, digits=3))", + # # ylim=(-1.1, 1.1), + # # xlim=(ymin, ymax)) - # plot!(ydisc, exact[i, :], - # label="Analytical", - # linewidth=2, - # linestyle=:dash) + # # plot!(ydisc, exact[i, :], + # # label="Analytical", + # # linewidth=2, + # # linestyle=:dash) - # end - # gif(anim, "pide_solution_mirrored_4pi.gif", fps=10) # Save the animation + # # end + # # gif(anim, "pide_solution_mirrored_4pi.gif", fps=10) # Save the animation end From 5ca4799e4de4dfeac39cf557420894a29c345a8d Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Sun, 3 Aug 2025 22:51:18 +0200 Subject: [PATCH 4/9] Reduce duplication of _euler_integral() for different dx types --- .../integral_expansion/integral_expansion.jl | 29 +++++++------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl index 1b8e777ed..7929f3e6a 100644 --- a/src/discretization/schemes/integral_expansion/integral_expansion.jl +++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl @@ -1,31 +1,24 @@ # use the trapezoid rule -function _euler_integral(II, s, jx, u, ufunc, dx::Number) #where {T,N,Wind,DX<:Number} +function _euler_integral(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} j, x = jx - if II[j] == 1 + if II[j] == 1 # recursively arrived at lower end of the domain return Num(0) end # unit index in direction of the derivative I1 = unitindex(ndims(u, s), j) # dx for multiplication Itap = [II - I1, II] - weights = [dx / 2, dx / 2] - return sym_dot(weights, ufunc(u, Itap, x)) + - _euler_integral(II - I1, s, jx, u, ufunc, dx) -end - -# Nonuniform dx -function _euler_integral(II, s, jx, u, ufunc, dx::AbstractVector) #where {T,N,Wind,DX<:Number} - j, x = jx - if II[j] == 1 - return Num(0) + if dx isa Number # Uniform dx + weights = [dx / 2, dx / 2] + elseif dx isa AbstractVector # Nonuniform dx + weights = fill(dx[II[j] - 1] / 2, 2) + else + error("Unsupported type of dx: $(type(dx))") end - # unit index in direction of the derivative - I1 = unitindex(ndims(u, s), j) - # dx for multiplication - Itap = [II - I1, II] - weights = fill(dx[II[j] - 1] / 2, 2) - + + # sym_do computes from II to II - I1, + # and recursive call computes from II - I1 to lower end of domain return sym_dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx) end From 46c8082d2647117f01c44267cdabdd97362398e0 Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Mon, 4 Aug 2025 09:45:00 +0200 Subject: [PATCH 5/9] Add transformation of [x...END]-integrals --- .../integral_expansion/integral_expansion.jl | 54 +++++++++++++++---- .../pde_system_transformation.jl | 16 +++--- 2 files changed, 54 insertions(+), 16 deletions(-) mode change 100644 => 100755 src/system_parsing/pde_system_transformation.jl diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl index 7929f3e6a..93f278d79 100644 --- a/src/discretization/schemes/integral_expansion/integral_expansion.jl +++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl @@ -22,11 +22,42 @@ function _euler_integral(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} return sym_dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx) end +function _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} + j, x = jx + if II[j] == length(s, x) # recursively arrived at upper end of the domain + return Num(0) + end + # unit index in direction of the derivative + I1 = unitindex(ndims(u, s), j) + # dx for multiplication + Itap = [II, II + I1] -function euler_integral(II, s, jx, u, ufunc) + if dx isa Number # Uniform dx + weights = [dx / 2, dx / 2] + elseif dx isa AbstractVector # Nonuniform dx + weights = fill(dx[II[j] + 1] / 2, 2) + else + error("Unsupported type of dx: $(type(dx))") + end + + # sym_do computes from II to II + I1, + # and recursive call computes from II + I1 to upper end of domain + return sym_dot(weights, ufunc(u, Itap, x)) + + _euler_integral_II_to_upper(II + I1, s, jx, u, ufunc, dx) +end + +# An integral from II to end of domain +function euler_integral(method::Symbol, II, s, jx, u, ufunc) j, x = jx dx = s.dxs[x] - return _euler_integral(II, s, jx, u, ufunc, dx) + if method == :lower_boundary_to_x + return _euler_integral(II, s, jx, u, ufunc, dx) + elseif method == :x_to_upper_boundary + return _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) + # error("Method :x_to_upper_boundary is not implemented for euler_integral.") + else + error("Unsupported method for euler_integral: $method") + end end # An integral across the whole domain (xmin .. xmax) @@ -47,13 +78,18 @@ end II::CartesianIndex, s::DiscreteSpace, depvars, indexmap, terms) ufunc(u, I, x) = s.discvars[u][I] - eulerrules = reduce(safe_vcat, - [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], - Num(x)))(u) => euler_integral( - Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc) - for x in ivs(u, s)] - for u in depvars], - init = []) + eulerrules = reduce(safe_vcat, [ + reduce(safe_vcat, [ + [ + # integrals from lower domain end to x: + Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(:lower_boundary_to_x, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), + # integrals from x to upper domain end: + Integral(x in DomainSets.ClosedInterval(Num(x), s.vars.intervals[x][2]))(u) => euler_integral(:x_to_upper_boundary, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), + # TODO: # any other arbitrary integrals??? + # TODO... + ] + for x in ivs(u, s)]) for u in depvars] ) + return eulerrules end diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl old mode 100644 new mode 100755 index 9c47960e5..78dbde0df --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -156,11 +156,13 @@ function descend_to_incompatible(term, v) end elseif op isa Integral if any(isequal(op.domain.variables), v.x̄) - euler = isequal( - op.domain.domain.left, v.intervals[op.domain.variables][1]) && - isequal(op.domain.domain.right, Num(op.domain.variables)) - whole = isequal( - op.domain.domain.left, v.intervals[op.domain.variables][1]) && + euler = ( # from lower domain boundary to x: + isequal(op.domain.domain.left, v.intervals[op.domain.variables][1]) && + isequal(op.domain.domain.right, Num(op.domain.variables))) || + (# from x to upper domain boundary: + isequal(op.domain.domain.left, Num(op.domain.variables)) && + isequal(op.domain.domain.right, v.intervals[op.domain.variables][2])) + whole = isequal(op.domain.domain.left, v.intervals[op.domain.variables][1]) && isequal(op.domain.domain.right, v.intervals[op.domain.variables][2]) if any([euler, whole]) u = arguments(term)[1] @@ -168,10 +170,10 @@ function descend_to_incompatible(term, v) @assert out==(nothing, false) "Integral $term must be purely of a variable, got $u. Try wrapping the integral argument with an auxiliary variable." return (nothing, nothing, false) else - throw(ArgumentError("Integration Domain only supported for integrals from start of iterval to the variable, got $(op.domain.domain) in $(term)")) + throw(ArgumentError("Integration Domain only supported for integrals across whole interval, or from an interval boundary (upper or lower) to the independent variable, got $(op.domain.domain) in $(term)")) end else - throw(ArgumentError("Integral must be with respect to the independent variable in its upper-bound, got $(op.domain.variables) in $(term)")) + throw(ArgumentError("Integral must be with respect to the independent variable in its upper- or lower-bound, got $(op.domain.variables) in $(term)")) end end for arg in SU.arguments(term) From eb8dd8cb55098f032e5094d0a9c95fb6c6cd042a Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Mon, 4 Aug 2025 10:11:07 +0200 Subject: [PATCH 6/9] Activate tests for [x..END]-integrals --- test/pde_systems/MOL_1D_Integration.jl | 130 ++++++++++++------------- 1 file changed, 64 insertions(+), 66 deletions(-) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index 1db4ecfdd..4f8d8b3c5 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -254,41 +254,40 @@ end # Discretize and compute numerical solution disc = MOLFiniteDifference([y => 50], t) - @test_broken (discretize(pde_system, disc) isa ODEProblem) - # prob = discretize(pde_system, disc) - # sol = solve(prob, Tsit5(), saveat=0.1) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) - # # Extract numerical solution at grid points - # ydisc = sol[y] - # tdisc = sol[t] - # usol = sol[u(t,y)] - - # # Calculate analytical solution at grid points - # # Since u(t,x) = x * (1 - exp(-t)) and y = -x, then u(t,y) = -y * (1 - exp(-t)) - # analytical_solution(t, y) = -y * (1 - exp(-t)) - # exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] - - # # Compare numerical and analytical solutions - # @test usol ≈ exact atol=0.1 - - # # # Create animated plot - # # using Plots # Add this import - # # anim = @animate for (i, t_val) in enumerate(tdisc) - # # plot(ydisc, usol[i, :], - # # label="Numerical", - # # linewidth=2, - # # xlabel="x", - # # ylabel="u(t,x)", - # # title="Solution at t = $(round(t_val, digits=3))", - # # ylim=(0, maximum(usol)*1.1), - # # xlim=(xmin, xmax)) + # Extract numerical solution at grid points + ydisc = sol[y] + tdisc = sol[t] + usol = sol[u(t,y)] + + # Calculate analytical solution at grid points + # Since u(t,x) = x * (1 - exp(-t)) and y = -x, then u(t,y) = -y * (1 - exp(-t)) + analytical_solution(t, y) = -y * (1 - exp(-t)) + exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.1 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(ydisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="x", + # ylabel="u(t,x)", + # title="Solution at t = $(round(t_val, digits=3))", + # ylim=(0, maximum(usol)*1.1), + # xlim=(xmin, xmax)) - # # plot!(xdisc, exact[i, :], - # # label="Analytical", - # # linewidth=2, - # # linestyle=:dash) - # # end - # # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation + # plot!(xdisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + # end + # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation end @testset "Test 04: PIDE with spatial derivative and source term - domain [0, 4π]" begin @@ -400,40 +399,39 @@ end # Discretize and compute numerical solution disc = MOLFiniteDifference([y => 200], t) - @test_broken (discretize(pde_system, disc) isa ODEProblem) - # prob = discretize(pde_system, disc) - # sol = solve(prob, Tsit5(), saveat=0.1) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) - # # Extract numerical solution at grid points - # ydisc = sol[y] - # tdisc = sol[t] - # usol = sol[u(t,y)] - - # # Calculate analytical solution at grid points - # # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) - # analytical_solution(t, y) = -sin(y) * cos(t) - # exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] - - # # Compare numerical and analytical solutions - # @test usol ≈ exact atol=0.36 - - # # # Create animated plot - # # using Plots # Add this import - # # anim = @animate for (i, t_val) in enumerate(tdisc) - # # plot(ydisc, usol[i, :], - # # label="Numerical", - # # linewidth=2, - # # xlabel="y", - # # ylabel="u(t,y)", - # # title="PIDE Solution [-4π,0] at t = $(round(t_val, digits=3))", - # # ylim=(-1.1, 1.1), - # # xlim=(ymin, ymax)) + # Extract numerical solution at grid points + ydisc = sol[y] + tdisc = sol[t] + usol = sol[u(t,y)] + + # Calculate analytical solution at grid points + # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) + analytical_solution(t, y) = -sin(y) * cos(t) + exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.36 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(ydisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="y", + # ylabel="u(t,y)", + # title="PIDE Solution [-4π,0] at t = $(round(t_val, digits=3))", + # ylim=(-1.1, 1.1), + # xlim=(ymin, ymax)) - # # plot!(ydisc, exact[i, :], - # # label="Analytical", - # # linewidth=2, - # # linestyle=:dash) + # plot!(ydisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) - # # end - # # gif(anim, "pide_solution_mirrored_4pi.gif", fps=10) # Save the animation + # end + # gif(anim, "pide_solution_mirrored_4pi.gif", fps=10) # Save the animation end From 0a9203f66aec44df48ba61566f528584c169f3bb Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Mon, 4 Aug 2025 11:08:31 +0200 Subject: [PATCH 7/9] Fix test 03b --- test/pde_systems/MOL_1D_Integration.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index 4f8d8b3c5..74a985b99 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -227,7 +227,7 @@ end @testset "Test 03b: Mirrored version of: Partial integro-differential equation with time derivative, with sys transformation" begin # Original equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t)) # Mirrored version with transformation y = -x, x = -y: - # Equation: ∂u/∂t + u(t,y) + ∫₀ˣ u(t,ξ) dξ = -x + x²/2 * (1 - exp(-t)) + # Equation: ∂u/∂t + u(t,y) + ∫_y ^⁰ u(t,ξ) dξ = -y + y²/2 * (1 - exp(-t)) # Initial condition: u(0,y) = 0 # Boundary condition: u(t,0) = 0 # Analytical solution: u(t,y) = -y * (1 - exp(-t)) @@ -242,8 +242,7 @@ end # Integral from y to 0 (note the reversed limits) Iy_to_0 = Integral(y in DomainSets.ClosedInterval(y, ymax)) - #eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t))] - eqs = [Dt(u(t,y)) + u(t,y) - Iy_to_0(u(t,y)) ~ -y + y^2/2*(1-exp(-t))] + eqs = [Dt(u(t,y)) + u(t,y) + Iy_to_0(u(t,y)) ~ -y + y^2/2*(1-exp(-t))] bcs = [ u(0,y) ~ 0.0, # Initial condition u(t,0) ~ 0.0] # Boundary condition (at y = 0, which corresponds to x = 0) @@ -276,13 +275,13 @@ end # plot(ydisc, usol[i, :], # label="Numerical", # linewidth=2, - # xlabel="x", - # ylabel="u(t,x)", + # xlabel="y", + # ylabel="u(t,y)", # title="Solution at t = $(round(t_val, digits=3))", # ylim=(0, maximum(usol)*1.1), - # xlim=(xmin, xmax)) + # xlim=(ymin, ymax)) - # plot!(xdisc, exact[i, :], + # plot!(ydisc, exact[i, :], # label="Analytical", # linewidth=2, # linestyle=:dash) From a7073b22454961f466531fef8e9385ee36b40dad Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Mon, 4 Aug 2025 11:28:13 +0200 Subject: [PATCH 8/9] Clarify naming --- .../generate_finite_difference_rules.jl | 27 ++++++++++++------- .../integral_expansion/integral_expansion.jl | 21 ++++++--------- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index 7ae8b4760..f4ce35b61 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -77,24 +77,29 @@ function generate_finite_difference_rules( # Spherical diffusion scheme spherical_diffusion_rules = generate_spherical_diffusion_rules( II, s, depvars, derivweights, bmap, indexmap, split_additive_terms(pde)) - integration_rules = vec(generate_euler_integration_rules( - II, s, depvars, indexmap, terms)) + + # Integration schemes for boundaries: [lower-bound, x] and [x, upper-bound] + integration_rules_from_boundary_to_x = generate_euler_integration_rules( + II, s, depvars, indexmap, terms) else central_deriv_rules_cartesian = [] advection_rules = [] nonlinlap_rules = [] spherical_diffusion_rules = [] mixed_deriv_rules_cartesian = [] - integration_rules = [] + integration_rules_from_boundary_to_x = [] end cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) - integration_rules = vcat(integration_rules, - vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) + # Integration schemes for boundaries: [lower-bound, upper-bound] + whole_domain_integration_rules = generate_whole_domain_integration_rules( + II, s, depvars, indexmap, terms) + return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), - vec(advection_rules), integration_rules) + vec(advection_rules), + vec(integration_rules_from_boundary_to_x), vec(whole_domain_integration_rules)) end function generate_finite_difference_rules( @@ -112,10 +117,12 @@ function generate_finite_difference_rules( advection_rules = [] nonlinlap_rules = [] spherical_diffusion_rules = [] - integration_rules = [] + integration_rules_from_boundary_to_x = [] + + whole_domain_integration_rules = generate_whole_domain_integration_rules( + II, s, depvars, indexmap, terms) - integration_rules = vcat(integration_rules, - vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), - vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) + vec(central_deriv_rules_cartesian), vec(advection_rules), + vec(integration_rules_from_boundary_to_x), vec(whole_domain_integration_rules)) end diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl index 93f278d79..3949f83e6 100644 --- a/src/discretization/schemes/integral_expansion/integral_expansion.jl +++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl @@ -1,5 +1,5 @@ # use the trapezoid rule -function _euler_integral(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} +function _euler_integral_lower_to_II(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} j, x = jx if II[j] == 1 # recursively arrived at lower end of the domain return Num(0) @@ -20,7 +20,7 @@ function _euler_integral(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} # sym_do computes from II to II - I1, # and recursive call computes from II - I1 to lower end of domain return sym_dot(weights, ufunc(u, Itap, x)) + - _euler_integral(II - I1, s, jx, u, ufunc, dx) + _euler_integral_lower_to_II(II - I1, s, jx, u, ufunc, dx) end function _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} j, x = jx @@ -51,7 +51,7 @@ function euler_integral(method::Symbol, II, s, jx, u, ufunc) j, x = jx dx = s.dxs[x] if method == :lower_boundary_to_x - return _euler_integral(II, s, jx, u, ufunc, dx) + return _euler_integral_lower_to_II(II, s, jx, u, ufunc, dx) elseif method == :x_to_upper_boundary return _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) # error("Method :x_to_upper_boundary is not implemented for euler_integral.") @@ -64,14 +64,11 @@ end function whole_domain_integral(II, s, jx, u, ufunc) j, x = jx dx = s.dxs[x] - if II[j] == length(s, x) - return _euler_integral(II, s, jx, u, ufunc, dx) - end - - dist2max = length(s, x) - II[j] + dist2max = length(s, x) - II[j] # distance to the end of the domain I1 = unitindex(ndims(u, s), j) Imax = II + dist2max * I1 - return _euler_integral(Imax, s, jx, u, ufunc, dx) + + return _euler_integral_lower_to_II(Imax, s, jx, u, ufunc, dx) end @inline function generate_euler_integration_rules( @@ -82,11 +79,9 @@ end reduce(safe_vcat, [ [ # integrals from lower domain end to x: - Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(:lower_boundary_to_x, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), + Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(:lower_boundary_to_x, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), # integrals from x to upper domain end: - Integral(x in DomainSets.ClosedInterval(Num(x), s.vars.intervals[x][2]))(u) => euler_integral(:x_to_upper_boundary, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), - # TODO: # any other arbitrary integrals??? - # TODO... + Integral(x in DomainSets.ClosedInterval(Num(x), s.vars.intervals[x][2]))(u) => euler_integral(:x_to_upper_boundary, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc) ] for x in ivs(u, s)]) for u in depvars] ) From 79cba5562571c12ea477bd0744932dcda41817c9 Mon Sep 17 00:00:00 2001 From: fabern <10245680+fabern@users.noreply.github.com> Date: Mon, 4 Aug 2025 12:01:55 +0200 Subject: [PATCH 9/9] Improve comments in tests --- test/pde_systems/MOL_1D_Integration.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index 74a985b99..b66217780 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -161,10 +161,11 @@ end end @testset "Test 03: Partial integro-differential equation with time derivative, with sys transformation" begin - # Equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t)) - # Initial condition: u(0,x) = 0 + # Equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t - τ)) with τ = 0.2 (alternatively τ = 0) + # On domain: x ∈ [0,2] and t ∈ [0,1] + # Initial condition: u(0,x) = x * (1 - exp(-τ)) # Boundary condition: u(t,0) = 0 - # Analytical solution: u(t,x) = x * (1 - exp(-t)) + # Analytical solution: u(t,x) = x * (1 - exp(-t - τ)) @parameters t, x @variables u(..) @@ -175,10 +176,8 @@ end Ix_0_to_x = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x - #eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t))] eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t - 0.2))] bcs = [ - # u(0,x) ~ 0.0, # Initial condition u(0,x) ~ x * (1 - exp(-0.2)), # Initial condition u(t,0) ~ 0.0] # Boundary condition domains = [ @@ -197,7 +196,7 @@ end usol = sol[u(t,x)] # Calculate analytical solution at grid points - analytical_solution(t, x) = x * (1 - exp(-(t+0.2)) ) + analytical_solution(t, x) = x * (1 - exp(-t - 0.2) ) exact = [analytical_solution(t_, x_) for t_ in tdisc, x_ in xdisc] # Compare numerical and analytical solutions @@ -225,9 +224,10 @@ end end @testset "Test 03b: Mirrored version of: Partial integro-differential equation with time derivative, with sys transformation" begin - # Original equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t)) + # Original equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t - τ)) with τ = 0 # Mirrored version with transformation y = -x, x = -y: # Equation: ∂u/∂t + u(t,y) + ∫_y ^⁰ u(t,ξ) dξ = -y + y²/2 * (1 - exp(-t)) + # On domain: y ∈ [-2, 0] and t ∈ [0,1] # Initial condition: u(0,y) = 0 # Boundary condition: u(t,0) = 0 # Analytical solution: u(t,y) = -y * (1 - exp(-t)) @@ -290,13 +290,13 @@ end end @testset "Test 04: PIDE with spatial derivative and source term - domain [0, 4π]" begin - # Equation: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) - # where h(x,t) = -sin(x)sin(t) + sin(x)cos(t) - cos(x)cos(t) + cos(t)(1 - cos(x)) + # Equation: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + # where h(x,t) = -sin(x)sin(t) + sin(x)cos(t) - cos(x)cos(t) + cos(t)(1 - cos(x)) # On domain: x ∈ [0,L]=[0,2π] and t ∈ [0,T]=[0,1] # Initial condition: u(0,x) = sin(x) # Boundary conditions: u(t,0) = 0, u(t,2π) = sin(2π)cos(t) = 0 # Analytical solution: u(t,x) = sin(x)cos(t) - + @parameters t, x @variables u(..) Dt = Differential(t)