PDEDiffusion1D

Diffusion 1D

import Pkg
Pkg.add("Plots")
using Plots

# Parameters
L = 1.0          # Length of domain
Nx = 100         # Number of spatial grid points
dx = L/(Nx-1)    # Spatial step
D  = 0.01        # Diffusion coefficient
dt = 0.0005      # Time step
Nt = 400         # Number of time steps
x = LinRange(0, L, Nx)

# Finite difference Laplacian (second derivative)
function laplacian(u, dx)
    du = similar(u)
    du[1] = 0.0           # Neumann BC: du/dx = 0 at boundaries
    du[end] = 0.0
    for i in 2:(length(u)-1)
        du[i] = (u[i-1] - 2u[i] + u[i+1]) / dx^2
    end
    return du
end

# Time stepping (Forward Euler)
function diffuse(u0, D, dx, dt, Nt)
    u = copy(u0)
    sol = [copy(u)]
    for n in 1:Nt
        u .= u .+ D*dt .* laplacian(u, dx)
        push!(sol, copy(u))
    end
    return sol
end

# Different initial conditions
u0_gaussian = exp.(-(x .- 0.5).^2 ./ 0.01)
u0_step     = x .< 0.5
u0_sine     = sin.(π*x)

# Solve
sol_gaussian = diffuse(u0_gaussian, D, dx, dt, Nt)
#sol_step     = diffuse(u0_step, D, dx, dt, Nt)
sol_sine     = diffuse(u0_sine, D, dx, dt, Nt)

# Plot snapshots
function plot_solution(sol, x; title_str="Diffusion")
    times = [1, round(Int,Nt/4), round(Int,Nt/2), Nt]
    plot(x, sol[1], label="t=0")
    for t in times[2:end]
        plot!(x, sol[t], label="t=$(round(t*dt, digits=2))")
    end
    xlabel!("x")
    ylabel!("u(x,t)")
    title!(title_str)
end

plot1 = plot_solution(sol_gaussian, x, title_str="Gaussian Initial Condition")
#plot2 = plot_solution(sol_step, x, title_str="Step Initial Condition")
plot3 = plot_solution(sol_sine, x, title_str="Sine Initial Condition")

plot(plot1, plot3, layout=(3,1), size=(700,900))

   Resolving package versions...
  No Changes to `/private/var/folders/p6/pwl670jd4sgf2p8z28gykgd00000gn/T/jl_0N9nw7/QuartoSandbox/Project.toml`
  No Changes to `/private/var/folders/p6/pwl670jd4sgf2p8z28gykgd00000gn/T/jl_0N9nw7/QuartoSandbox/Manifest.toml`