Визуализация решения дифференциального уравнения с анимацией в Юлии - PullRequest
Новая
       30

Визуализация решения дифференциального уравнения с анимацией в Юлии

1 голос
/ 20 апреля 2019

Я решил систему двух дифференциальных уравнений второго порядка, используя реализацию метода Эйлера в Джулии. Приведенный ниже код показывает, как был вызван метод Эйлера для решения рассматриваемой системы. 

θ1 = 1.1900518004210798; θ2 = 0.3807445263738167
f(t, y) = [y[2], -2(y[1] - θ1) - 4y[2] + 0.5sin(3pi*t),
           y[4], -2(y[3] - θ2) - 4(y[4] + abs(y[2])) + 0.5sin(3pi*t)]
y0 = [pi/2, 0, pi/6, 0]; t0 = 0; tFinal = 50; h = 0.001
res = euler(f, y0, t0, tFinal, h)

Результат, res, представляет собой вектор из четырех чисел 

1.18798735437173    
-0.0458294959470722  
0.31530569612003573 
-0.049213402534541074

Первая запись - это угол, который образует нижний сегмент линии с осью X, а третья запись - это угол, который образуют два отрезка линии друг с другом (см. Рисунок ниже).

enter image description here

Для создания этого сюжета я назвал plot_robotarm([res[1], res[3]]), который реализован в соответствии с приведенным ниже кодом. 

function plot_robotarm(thetav)
    # Plots a robotarm with angles thetav
    R = 1;
    xv=zeros(length(thetav)+1)
    yv=zeros(length(thetav)+1)
    for i in 1:length(thetav)
        xv[i+1]=xv[i]+R*cos(thetav[i])
        yv[i+1]=yv[i]+R*sin(thetav[i])
    end
    # Plot with colors 
    opts = (:circle, 10, 1., :blue, stroke(7, 1., :red))
    plt = plot(xv, yv,
               marker = opts,
               c = :red,
               w = 5,
               legend = false,
               xlims = (0, 2.0),
               ylims = (0, 2.0))
    display(plt)
end

Как я могу создать анимацию, которая визуализирует, как последовательные итерации метода Эйлера заставляют манипулятор робота (то есть два отрезка) двигаться к конечной точке в t = 50? (Мне не нужно составлять график каждой итерации, достаточно того, чтобы она создавала анимацию, которая фиксирует движение руки робота.)

1 Ответ

2 голосов
/ 21 апреля 2019

Вы можете использовать анимационные функции ffmpeg и Luxor.jl для создания анимированного GIF.Функцию фрейма необходимо изменить, чтобы отразить графическое отображение каждого шага в вашей программе.Смотрите документы для Луксора больше. 

using Luxor
using Colors
using BoundaryValueDiffEq

# constants for differential equations and movie
const g = 9.81
const L = 1.0                         # pendulum length in meters
const bobd = 0.10                     # pendulum bob diameter in meters
const framerate = 50.0                # intended frame rate/sec
const t0 = 0.0                        # start time (s)
const tf = 2.3                        # end simulation time (s)
const dtframe = 1.0/framerate         # time increment per frame
const tspan = LinRange(t0, tf, Int(floor(tf*framerate)))  # array of time points in animation

const bgcolor = "black"               # gif background
const leaderhue = (0.80, 0.70, 0.20)  # gif swing arm hue light gold
const hslcolors = [HSL(col) for col in (distinguishable_colors(
                   Int(floor(tf*framerate)+3),[RGB(1,1,1)])[2:end])]
const giffilename = "pendulum.gif"    # output file

# differential equations copied from docs of DifferentialEquations.jl
simplependulum!(du, u, p, t) = (θ=u[1]; dθ=u[2]; du[1]=dθ; du[2]=-(g/L)*sin(θ))
bc1!(residual, u, p, t) = (residual[1] = u[div(end,2)][1] + pi/2; residual[2] = u[end][1] - pi/2)
bvp1 = TwoPointBVProblem(simplependulum!, bc1!, [pi/2,pi/2], (tspan[1],tspan[end]))
sol2 = solve(bvp1, GeneralMIRK4(), dt=dtframe) 

# movie making background
backdrop(scene, framenumber) = background(bgcolor)

function frame(scene, framenumber)
    u1, u2 = sol2.u[framenumber]
    y, x = L*cos(u1), L*sin(u1)
    sethue(leaderhue)
    poly([Point(-4.0, 0.0), Point(4.0, 0.0),
          Point(160.0x,160.0y)], :fill)
    sethue(Colors.HSV(framenumber*4.0, 1, 1))
    circle(Point(160.0x,160.0y), 160bobd, :fill)
    text(string("frame $framenumber of $(scene.framerange.stop)"),
        Point(0.0, -190.0),
        halign=:center)
end

muv = Movie(400, 400, "Pendulum Demo", 1:length(tspan))
animate(muv, [Scene(muv, backdrop),
              Scene(muv, frame, easingfunction=easeinoutcubic)],
              creategif=true, pathname=giffilename)
...