Skip to content

Commit 97b28fe

Browse files
Create ODE_Nonlinear_system_of_reactions_with_an_analytical_solution
1 parent aead926 commit 97b28fe

File tree

1 file changed

+28
-0
lines changed

1 file changed

+28
-0
lines changed
Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
using DifferentialEquations
2+
using Plots
3+
using SpecialFunctions
4+
5+
function nonLinChem(dy,y,p,t)
6+
dy[1] = -y[1]
7+
dy[2] = y[1]-(y[2])^2
8+
dy[3] = (y[2])^2
9+
end
10+
y0 = [1.0;0.0;0.0]
11+
tspan = (0.0,20)
12+
prob = ODEProblem(nonLinChem,y0,tspan)
13+
sol = solve(prob)
14+
15+
Analytical(t) = [exp(-t);
16+
(2sqrt(exp(-t))besselk(1,2sqrt(exp(-t)))-2besselk(1,2)/besseli(1,2)*sqrt(exp(-t))besseli(1,2sqrt(exp(-t))))/(2besselk(0,2sqrt(exp(-t)))+(2besselk(1,2)/besseli(1,2))besseli(0,2sqrt(exp(-t))));
17+
-exp(-t)+1+(-2sqrt(exp(-t))*besselk(1,2sqrt(exp(-t)))+sqrt(exp(-t))*besseli(1,2sqrt(exp(-t)))*2besselk(1,2)/besseli(1,2))/(2besselk(0,2sqrt(exp(-t)))+2besselk(1,2)/besseli(1,2)*besseli(0,2sqrt(exp(-t))))]
18+
19+
t=0:.1:20
20+
p = 1
21+
plot(sol,vars=(0,p))
22+
plot!(t,(x->Analytical(x)[p]).(t))
23+
p = 2
24+
plot!(sol,vars=(0,p))
25+
plot!(t,(x->Analytical(x)[p]).(t))
26+
p = 3
27+
plot!(sol,vars=(0,p))
28+
plot!(t,(x->Analytical(x)[p]).(t))

0 commit comments

Comments
 (0)