# 有关ODE求解的问题

#1

module msh

using PyPlot

export plotMesh

mutable struct Mesh
coords

``````function Mesh(n1,n2)
@assert (n1 > 0 && n2 > 0) "Mesh arguments must be positive"
coords = createCoords(n1,n2)
new(coords)
end
``````

end

function createCoords(n1,n2)

#Fill coords array
coords = zeros(2,n1n2);
cnt = 1
for i = 1:n2
for j = 1:n1
coords[1:2,cnt] = [j-1,i-1]
cnt = cnt+1
end
end
r = 1/(4
n1)*rand(size(coords)[1], size(coords)[2])
coords = coords + r #“random” pertubation of grid
return coords
end

function plotMesh(m::Mesh,c::String)

``````#Plot coordinates
for i = 1:size(m.coords)[2]
plot(m.coords[1,i],m.coords[2,i],marker="o",color =c, markersize = 10)
end
``````

end
end

using PyPlot
using ODE
using LinearAlgebra
###########################################################################################

# INPUT

###########################################################################################

n = 10 #number of particles will be nn
N = n
n
#set up mesh
m = msh.Mesh(n,n)
npull = n #number of elements being pulled at a constant speed
nhold = n #number of elements being held at their initial position

#Set time interval [0,maxT]
maxT = 4

#Set up matrix of spring constants
global K
K = zeros(N,N)
for i = 1:N
for j = i+1:N
K[i,j] = 2^(-norm(m.coords[:,i] - m.coords[:,j],2)+1) #compute 2-norm because why not
K[j,i] = K[i,j]
end
end
#strengthen bond to the boundary nodes
for i = 1:nhold
K[i,:] *= 10
K[:,i] *= 10
end
for i = N-npull+1:N
K[i,:] *= 10
K[:,i] *= 10
end

##################################################################################################

#count the number of bonds each element has
global nbonds = Dict{Float64, Array{Float64,2}}()
nbonds[0.] = N*ones(N,1)

c1 = [4 2] #choose speed in x and y direction
function f(t,x)
global K
H = ones(size(K)[1],size(K)[2])
xh = [m.coords[1,1:nhold]’ m.coords[2,1:nhold]’;
hcat(x[1:N-npull-nhold], x[N-npull-nhold+1:2*(N-npull-nhold)]);
(c1[1]*t+m.coords[1,N-npull+1:N])’ (c1[2]*t+m.coords[2,N-npull+1:N])’]’
for i = 1:N
for j = i+1:N
if norm(xh[:,i]-xh[:,j]) > 10
H[i,j] = H[j,i] = 0
end
#make particles repulse each other/try to enforce order
#confuses the ode solver for some reason
#if norm(xh[:,i] - xh[:,j]) < 0.1 && H[i,j] > 0
# H[i,j] = -H[i,j]
# H[j,i] = H[i,j]
#end
end
end
K = K.*H
#recalculate the number of bonds
global nbonds
nbonds[t] = sum(K.!=0, 2)

#Create stiffness matrix
A = zeros(N-npull-nhold,N-npull-nhold)
for i = 1:N-npull-nhold
A[i,i] = -sum(K[i+nhold,1:N]) #diagonal is special
for j = 1:N-npull-nhold
if(i!=j) A[i,j] = K[i+nhold,j+nhold] end
end
end
#Create “forcing” term
B = zeros(2*(N-npull-nhold),1)
for i = 1:N-npull-nhold
B[i] = sum((c1[1]*t + m.coords[1,N-npull+1:N]).K[i+nhold,N-npull+1:N]) + sum(K[i+nhold, 1:nhold].m.coords[1,1:nhold])
B[i+N-npull-nhold] = sum((c1[2]t + m.coords[2,N-npull+1:N]).K[i+nhold,N-npull+1:N]) + sum(K[i+nhold, 1:nhold].m.coords[2,1:nhold])
end
y = [A
x[1:N-npull-nhold]+B[1:N-npull-nhold]; A
x[N-npull-nhold+1:2
(N-npull-nhold)]+B[N-npull-nhold+1:2
(N-npull-nhold)]]’’
return [x[2
(N-npull-nhold)+1:end]; y]
end

#Solve system of second order ODE using ode45
y0 = [m.coords[1,nhold+1:N-npull]; m.coords[2,nhold+1:N-npull]; zeros(2*(N-npull-nhold))] #initial conditions
#initial positions given by mesh, initial speed set to zero
tout,yout = ode45(f, y0,[0:0.1:maxT;])
ys = hcat(yout…)’
@printf “Number of timesteps: %d\n” size(tout)[1]

#Plot position of particles at time step
cnt = 1
minh = minimum(tout[2:end] - tout[1:end-1])
diff = maximum(tout[2:end] - tout[1:end-1])/minh
if diff > 5 i_step = floor(diff/5) end
for tstep = 1:i_step:size(tout)[1] - 1
currenth = tout[tstep+1] - tout[tstep]
for j =1:i_step:ceil(currenth/minh)
broken_t = (nbonds[tout[tstep]] .< 50)[:,1]
yh = [m.coords[1,1:nhold]’ m.coords[2,1:nhold]’;
vcat(ys[tstep, 1:N-npull-nhold], ys[tstep, N-npull-nhold+1:2*(N-npull-nhold)])’;
(c1[1]*tout[tstep]+m.coords[1,N-npull+1:N])’ (c1[2]*tout[tstep]+m.coords[2,N-npull+1:N])’]
plot([-5,20],[-5,25], markersize = 0, linestyle = “None”) #sloppy way of setting axis sizes
plot(yh[ broken_t,1], yh[ broken_t,2], marker=“o”, markersize=10, color = “r”, linestyle = “None”)
plot(yh[~broken_t,1], yh[~broken_t,2], marker=“o”, markersize=10, color = “b”, linestyle = “None”)
savefig(string(“Plots/plot”, cnt))
clf()
cnt = cnt+1
end
end

ERROR: LoadError: MethodError: no method matching +(::Float64, ::Array{Float64,1})
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any…) at operators.jl:502
+(::Float64, ::Float64) at float.jl:395
+(::PyCall.PyObject, ::Any) at C:\Users\好兔子.julia\packages\PyCall\rUul9\src\pyoperators.jl:14

Stacktrace:
[1] f(::Float64, ::Array{Float64,1}) at D:\新建文件夹\Peridynamics-master\test.jl:91
[2] hinit(::typeof(f), ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Float64, ::Float64) at
C:\Users\好兔子.julia\packages\ODE\AfTL\src\ODE.jl:92
[3] #oderk_adapt#13(::Float64, ::Float64, ::Function, ::Float64, ::Float64, ::Int64, ::Symbol, ::Fu
nction, ::typeof(f), ::Array{Float64,1}, ::Array{Float64,1}, ::ODE.TableauRKExplicit{:dopri,7,Ration
al{Int64}}) at C:\Users\好兔子.julia\packages\ODE\AfTL\src\runge_kutta.jl:283
7,Rational{Int64}}) at C:\Users\好兔子.julia\packages\ODE\AfTL\src\runge_kutta.jl:244
[5] #ode45_dp#7 at C:\Users\好兔子.julia\packages\ODE\AfTL\src\runge_kutta.jl:217 [inlined]
[6] ode45_dp at C:\Users\好兔子.julia\packages\ODE\AfTL\src\runge_kutta.jl:217 [inlined]
[7] #ode45#8 at C:\Users\好兔子.julia\packages\ODE\AfTL\src\runge_kutta.jl:219 [inlined]
[8] ode45(::Function, ::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\好兔子.julia\packages\OD
E\AfTL\src\runge_kutta.jl:219
[9] top-level scope at none:0
[11] (::getfield(Atom, Symbol("##128#133")){String,String,Module})() at C:\Users\好兔子.julia\pack
ages\Atom\Fha1N\src\eval.jl:120
[12] withpath(::getfield(Atom, Symbol("##128#133")){String,String,Module}, ::String) at C:\Users\好

[13] withpath at C:\Users\好兔子.julia\packages\Atom\Fha1N\src\eval.jl:46 [inlined]
[14] #127 at C:\Users\好兔子.julia\packages\Atom\Fha1N\src\eval.jl:117 [inlined]
[15] hideprompt(::getfield(Atom, Symbol("##127#132")){String,String,Module}) at C:\Users\好兔子.ju
lia\packages\Atom\Fha1N\src\repl.jl:84
[16] macro expansion at C:\Users\好兔子.julia\packages\Atom\Fha1N\src\eval.jl:116 [inlined]
in expression starting at D:\新建文件夹\Peridynamics-master\test.jl:134

#2

#3

#4

# Current status of the ODE.jl project

This project is deprecated in favor of DifferentialEquations.jl and its ODE solvers OrdinaryDiffEq.jl. This library is in “maitanance mode”, meaning that it is being upgraded with each Julia version, but not seeing active feature development. ODE.jl contains the basic functionality that was moved here when the package was originally moved from Base. Although quite poorly tested, at least some of the functionality is quite reliable. Use at your own risk.

#5

#6

#7