Commit 4d644a43 authored by William Wright's avatar William Wright
Browse files

Merge branch 'main' of https://git.physics.byu.edu/justinct/Final-Project into main

parents 2d06cca0 65804372
......@@ -18,18 +18,35 @@ function lorenz!(du, u, p, t)
end
function compare(r1, v1, r2, v2)
x1 = r1[:][1]
y1 = r1[:][2]
z1 = r1[:][3]
xp1 = v1[:][1]
yp1 = v1[:][2]
zp1 = v1[:][3]
x2 = r2[:][1]
y2 = r2[:][2]
z2 = r2[:][3]
xp2 = v2[:][1]
yp2 = v2[:][2]
zp2 = v2[:][3]
x1 = []
y1 = []
z1 = []
x2 = []
y2 = []
z2 = []
xp1 = []
yp1 = []
zp1 = []
xp2 = []
yp2 = []
zp2 = []
for i in 1:100
push!(x1, r1[i][1])
push!(y1, r1[i][2])
push!(z1, r1[i][3])
push!(x2, r2[i][1])
push!(y2, r2[i][2])
push!(z2, r2[i][3])
push!(xp1, v1[i][1])
push!(yp1, v1[i][2])
push!(zp1, v1[i][3])
push!(xp2, v2[i][1])
push!(yp2, v2[i][2])
push!(zp2, v2[i][3])
end
xres = x2 .- x1
yres = y2 .- y1
......@@ -77,17 +94,28 @@ function timelorenzrun(tend, algorithm, initconds, params)
end
function getR(sol)
x = sol[1][:]
y = sol[2][:]
z = sol[3][:]
r = [x, y, z]
r = []
for i in 1:100
push!(r, sol(i))
end
return r
end
function getV(sol)
x = sol[1][:]
y = sol[2][:]
z = sol[3][:]
r = []
for i in 1:100
push!(r, sol(i))
end
x = []
y = []
z = []
for i in 1:100
push!(x, r[i][1])
push!(y, r[i][2])
push!(z, r[i][3])
end
xp = params[1].*(y .- x)
yp = x.*(params[2] .- z) .- y
......@@ -105,8 +133,16 @@ sol4 = timelorenzrun(tend, QNDF(), u0, params);
sol5 = timelorenzrun(tend, ImplicitDeuflhardExtrapolation(max_order=7,min_order=4,init_order=4,sequence=:bulirsch), u0, params);
sol6 = timelorenzrun(tend, ExtrapolationMidpointDeuflhard(max_order=7,min_order=4,init_order=4,sequence=:bulirsch,threading=false), u0, params);
convert(Array, sol1)
convert(Array, sol2)
convert(Array, sol3)
convert(Array, sol4)
convert(Array, sol5)
convert(Array, sol6)
solList = [sol1, sol2, sol3, sol4, sol5, sol6]
# Getting position and velocity data from solutions
Rs = []
Vs = []
......@@ -131,13 +167,13 @@ function plotRes(res)
yp = res[2][2]
zp = res[2][3]
plotx = plot(x, title= "Residual in X over time")
ploty = plot(y, title= "Residual in Y over time")
plotz = plot(z, title= "Residual in Z over time")
plotxp = plot(xp, title= "Residual in X velocity over time")
plotyp = plot(yp, title= "Residual in Y velocity over time")
plotzp = plot(zp, title= "Residual in Z velocity over time")
plotx = plot(x, title= "Residual in X ")
ploty = plot(y, title= "Residual in Y")
plotz = plot(z, title= "Residual in Z")
plotxp = plot(xp, title= "Residual in X velocity")
plotyp = plot(yp, title= "Residual in Y velocity")
plotzp = plot(zp, title= "Residual in Z velocity")
plot(plotx, ploty, plotz, plotxp, plotyp, plotzp, layout = (3,2))
plot(plotx, ploty, plotz, plotxp, plotyp, plotzp, layout = (3,2), xlabel = "Time Step", ylabel = "Residual")
end
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment