Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

自分用:Newmark-β法

More than 5 years have passed since last update.

早速コード

function newmarkβ(input, Δt, ω, h, β=0.25, ε=0.000001)
  output = zeros(input)

  u̇ₓ = u̇₁ = zeros(input[:, 1])
  uₓ = u₁ = output[:, 1] = zeros(input[:, 1])
  üₓ = ü₁ = -input[:, 1] - 2 * h * ω * u̇₁ - ω ^ 2 * u₁

  for x=1:length(input)-1
    ǖₓ₊₁ = üₓ
    while true
      u̇ₓ₊₁ = u̇ₓ + 0.5 * (üₓ + ǖₓ₊₁) * Δt
      uₓ₊₁ = uₓ + u̇ₓ * Δt + (0.5 - β) * üₓ * Δt ^ 2 + β * ǖₓ₊₁ * Δt ^ 2

      üₓ₊₁ = -input[:, x+1] - 2 * h * ω * u̇ₓ₊₁ - ω ^ 2 * uₓ₊₁

      if norm(ǖₓ₊₁ - üₓ₊₁) < ε
        üₓ = üₓ₊₁
        u̇ₓ = u̇ₓ₊₁
        uₓ = uₓ₊₁
        output[:, x+1] = uₓ₊₁
        break
      end
      ǖₓ₊₁ = üₓ₊₁
    end
  end
  return output
end

t = transpose(linspace(0.0, 100, 1000))
Δt = t[2]
input = sin(t)
res=newmarkβ(input, Δt, 1.0, 0.1)
using Gadfly
plot(x=t, y=res)

擬似コードに見せかけたJuliaで実行可能なコード
ユニコードプログラミングをネタ的にやってみた
普通に変位応答を戻り値で返しているが、実際にはobserverを受け取る実装にしたほうが良いだろう。
ここでは線形の範囲(多自由度にも一応対応しているはず)で解いているので全くNewmark-βの意味が無い。
非線形応答をやるなら剛性と減衰を関数か何かで渡す必要がある。コルーチンがあれば楽かな多分

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away