有关ODE求解的问题

haotuzi · 2018 年10 月 5 日 12:25

初次学习Julia，想把别人的老版本的近场动力代码在1.0上运行。。。。

代码
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/(4n1)*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 = nn
#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 = [Ax[1:N-npull-nhold]+B[1:N-npull-nhold]; Ax[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
[4] oderk_adapt(::Function, ::Array{Float64,1}, ::Array{Float64,1}, ::ODE.TableauRKExplicit{:dopri,
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
[10] include_string(::Module, ::String, ::String) at .\loading.jl:1002
[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\好
兔子.julia\packages\CodeTools\hB4Hy\src\utils.jl:30
[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]
[17] (::getfield(Atom, Symbol(“##126#131”)){Dict{String,Any}})() at .\task.jl:85
in expression starting at D:\新建文件夹\Peridynamics-master\test.jl:134

Gnimuc · 2018 年10 月 5 日 13:30

学会看 Stacktrace，这行的错误是：

需要将那行的 + 改成 .+，即 element-wise 操作（老版本要求不严格，1.0 要求显式写 .）。

Scheme · 2018 年10 月 6 日 04:20

不要用ODE.jl。

woclass · 2018 年10 月 6 日 08:35

补充

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.

ref: JuliaDiffEq/ODE.jl: Assorted basic Ordinary Differential Equation solvers

所以该改用以下包

JuliaDiffEq/DifferentialEquations.jl: Julia suite for high-performance solvers of differential equations

This is a suite for numerically solving differential equations in Julia. The purpose of this package is to supply efficient Julia implementations of solvers for various differential equations

haotuzi · 2018 年10 月 6 日 09:23

厄，那用哪个呢？

woclass · 2018 年10 月 6 日 10:11

用这个，原因见楼上

JuliaDiffEq/DifferentialEquations.jl

haotuzi · 2018 年10 月 6 日 10:16

非常感谢，我试试