JuliaScripts/ode_solution.jl at main · shenge86/JuliaScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
# ============================================================
#  Solving dy/dt = -1000y + 3000 - 2000*exp(-t),  y(0) = 0
#
#  Analytical solution:
#    y(t) = 3 - (2000/999)*exp(-t) - (997/999)*exp(-1000t)
#
#  We also solve numerically with DifferentialEquations.jl
#  using a stiff-aware solver (Rodas5), then compare both.
# ============================================================

using DifferentialEquations
using Plots

# ----------------------------------------------------------
# 1. Define the ODE right-hand side
# ----------------------------------------------------------
function ode!(du, u, p, t)
    du[1] = -1000*u[1] + 3000 - 2000*exp(-t)
end

# ----------------------------------------------------------
# 2. Set up and solve numerically
#    Rodas5 is an excellent stiff solver (5th-order Rosenbrock)
# ----------------------------------------------------------
u0    = [0.0]           # initial condition y(0) = 0
tspan = (0.0, 0.01)     # transient dies out very fast (rates 1 and 1000)
prob  = ODEProblem(ode!, u0, tspan)
sol   = solve(prob, Rodas5(), reltol=1e-10, abstol=1e-12)

# ----------------------------------------------------------
# 3. Analytical solution
# ----------------------------------------------------------
y_exact(t) = 3.0 - (2000/999)*exp(-t) - (997/999)*exp(-1000*t)

# Dense time grid for a smooth analytical curve
t_dense  = range(0.0, 0.01, length=1000)
y_dense  = y_exact.(t_dense)

# ----------------------------------------------------------
# 4. Plot
# ----------------------------------------------------------
p1 = plot(
    t_dense, y_dense,
    label  = "Analytical y(t)",
    lw     = 2.5,
    color  = :royalblue,
    xlabel = "t",
    ylabel = "y(t)",
    title  = "dy/dt = -1000y + 3000 - 2000·exp(-t),   y(0) = 0",
    legend = :bottomright,
    grid   = true,
    framestyle = :box,
)

plot!(p1,
    sol.t, first.(sol.u),
    label      = "Numerical (Rodas5)",
    lw         = 1.5,
    linestyle  = :dash,
    color      = :crimson,
    markershape = :circle,
    markersize  = 3,
)

hline!(p1, [3.0],
    label     = "Steady state y = 3",
    lw        = 1.2,
    linestyle = :dot,
    color     = :gray40,
)

# ----------------------------------------------------------
# 5. Zoom-in: show the ultra-fast initial transient
# ----------------------------------------------------------
t_zoom  = range(0.0, 0.002, length=500)
y_zoom  = y_exact.(t_zoom)

p2 = plot(
    t_zoom, y_zoom,
    label  = "Analytical (zoomed)",
    lw     = 2.5,
    color  = :royalblue,
    xlabel = "t",
    ylabel = "y(t)",
    title  = "Initial transient (t ∈ [0, 0.002])",
    legend = :bottomright,
    grid   = true,
    framestyle = :box,
)

hline!(p2, [3.0],
    label     = "Steady state y = 3",
    lw        = 1.2,
    linestyle = :dot,
    color     = :gray40,
)

# ----------------------------------------------------------
# 6. Combine and save
# ----------------------------------------------------------
fig = plot(p1, p2, layout = (2, 1), size = (800, 700), dpi = 150)
savefig(fig, "ode_solution.png")

println("Plot saved to ode_solution.png")
println()
println("Verification — exact solution at selected points:")
println("  y(0)     = ", round(y_exact(0.0),     digits=8), "  (should be 0)")
println("  y(0.001) = ", round(y_exact(0.001),   digits=8))
println("  y(0.005) = ", round(y_exact(0.005),   digits=8))
println("  y(0.01)  = ", round(y_exact(0.01),    digits=8))
println("  y(∞) ≈   ", round(y_exact(1.0),      digits=8), "  (→ 3)")