Test-2021-5-26-no1

Test-2021-5-26-no1-orginal.ipynb

Test-2021-5-26-no1-revised.ipynb

original から revised への変更点

もとの話題はツイッターより: https://twitter.com/genkuroki/status/1397762508010651654

--- Test-2021-5-26-no1-orginal.jl	2021-05-27 12:36:58.800973600 +0900
+++ Test-2021-5-26-no1-revised.jl	2021-05-27 12:37:16.652891900 +0900
@@ -17,6 +17,7 @@
 # +
 #May 26th, 2021
 #n重振り子を解きたいJulia
+#Revised: May 27th, 2021
 
 using LinearAlgebra,DifferentialEquations,Plots,StaticArrays
 
@@ -25,60 +26,41 @@
 g  = 9.8   # m/s²
 Δl = 0.0   # m
 Δm = 0.0   # kg
-p = (;n, g, Δl, Δm)
+# See https://github.com/JuliaArrays/StaticArrays.jl
+A = zeros(MMatrix{n,n})
+B = zeros(MMatrix{n,n})
+U = SVector{n}(-ones(n))
+L = SVector{n}(@. 1.0 + (1:n)*Δl)
+M = SVector{n}(@. 1.0 + (1:n)*Δm)
+p = (; n, g, Δl, Δm, A, B, U, L, M)
 
-function Npendulum(du, u, p, t)
-    n,g,Δl,Δm = p
-    L = []
-    M = []
-    U = -ones(n)
-    
-#setting of L,M
-for ii=1:n
-    lᵢ = 1.0 + ii*Δl
-    mᵢ = 1.0 + ii*Δm
-    push!(L, lᵢ)
-    push!(M, mᵢ)
-end
+function Npendulum!(du, u, p, t)
+    n, g, Δl, Δm, A, B, U, L, M = p
 
 #係数行列Aの初期化
-A = zeros(n,n)
     #対角成分
         for i=1:n
-            Mᵢᵢ = 0.0
-            for j=i:n
-                Mᵢᵢ += M[j]
-            end
+            Mᵢᵢ = sum(M[j] for j in i:n)
             A[i,i] = L[i]*L[i]*Mᵢᵢ
         end    
     #非対角成分
         for i=1:n
             for j=(i+1):n
-                Mᵢⱼ = 0.0
-                for k=j:n
-                    Mᵢⱼ += M[k]
-                end
+                Mᵢⱼ = sum(M[k] for k=j:n)
             A[i,j] = L[i]*L[j]*cos(u[i]-u[j])*Mᵢⱼ
             A[j,i] = A[i,j]
             end
         end    
 #係数行列Bの初期化
-B = zeros(n,n)
     #対角成分
         for i=1:n
-            Mᵢᵢ = 0.0
-            for j=i:n
-                Mᵢᵢ += M[j]
-            end
+            Mᵢᵢ = sum(M[j] for j=i:n)
             B[i,i] = g*L[i]*sin(u[i])*Mᵢᵢ
         end
     #非対角成分
         for i=1:n
             for j=(i+1):n
-                Mᵢⱼ = 0.0
-                for k=j:n
-                    Mᵢⱼ += M[k]
-                end
+                Mᵢⱼ = sum(M[k] for k=j:n)
             B[i,j] = L[i]*L[j]+u[n+j]^2*sin(u[i]-u[j])*Mᵢⱼ
             B[j,i] = L[j]*L[i]+u[n+i]^2*sin(u[j]-u[i])*Mᵢⱼ
             end
@@ -106,7 +88,7 @@
     u₀ = [θ ; ω]
     tinterval = (0.0,Lx)
 
-    prob2 = ODEProblem(Npendulum,u₀,tinterval,p)
+    prob2 = ODEProblem(Npendulum!,u₀,tinterval,p)
     @time sol = solve(prob2, Tsit5(), reltol=1e-8, abstol=1e-8)
     @time sol = solve(prob2, Tsit5(), reltol=1e-8, abstol=1e-8)
     @time sol = solve(prob2, Tsit5(), reltol=1e-8, abstol=1e-8)
@@ -130,4 +112,3 @@
 plot(Plt...; size=(xlength,ylength))
 # -
 -