Made a simple damped harmonic oscillator but need to make it driven instead of simple?

I’ve been writing a simple 3d spring for the last week or so now, and i’ve managed to get it to work when the target is being updated (however, for simple springs, its not supposed to work when the target is updated, or if it does, it won’t work how it should, that’s why it needs to be driven) however while the target is moving then stops the velocity of the spring suddenly increases, and if the target goes fast the slow, the spring doesn’t keep any of its momentum?

Also, I’m pretty sure i am using a few too many variables for my spring and i’m wondering if theres a way to improve the performance of it in the process?

This is my code, so far:

wait(0.5)

local target = workspace.Target
local current = workspace.Current

local function create(c, A, B)
	local s = {}
	
	--A for decay rate, B for frequency
	s.A = A
	s.B = B
	
	--s is the starting value, t is the target value, and c is the current value
	s.s = c
	s.t = c
	s.c = c
	
    --x is time elapsed since target not changed, v is velocity, p is previous position
	s.x = 0
	s.v = Vector3.new()
	s.p = Vector3.new()
	
	s.update = function(x)
        --set the total time elapsed
		s.x = s.x + x
        --set the velocity based on change of target position
		s.v = (s.t - s.p)
		s.p = s.t
        --the function which updates the spring
		s.c = (s.s - s.t) * math.cos(s.B * s.x) * 2^(-s.A * s.x) + s.t
        --when the target is moving
		if s.v ~= Vector3.new() then
			s.x = 0
			s.s = s.c
		end
	end
	return s
end

--creating the spring
local spring = create(Vector3.new(0, 0, 1.05), 2, 3)
current.Position = spring.c

--updating the spring
game:GetService("RunService").RenderStepped:Connect(function(s)
	spring.t = target.Position
	spring.update(s)
	current.Position = spring.c
end)

I’m pretty sure my method of updating the spring is also having a massive performance hit on the gpu, so making sure this is efficient as possible is my aim.

And before you think im talking about physical actual springs/constraints i’m not, i’m talking about similar spring mechanics found in fps games, such as phantom forces.

Here is a video to explain what I mean.

3 Likes

The problem is that you’re using the equation for cases with no external forces present. Since you’re dragging around the target in various ways, in real life that would induce a force on the spring which your equation doesn’t account for. Right now you are essentially using a model for a damped harmonic oscillator, whereas you need to model a driven harmonic oscillator.

I would look into using PID controllers. While they may seem somewhat arcane, they’re doing essentially the same thing, and in my opinion they’re much simpler to use than complicated differential equations (what a concept, I know). Plus, they perform much faster than anything involving trig functions and exponentials.

4 Likes

Are there any resources that look at code based driven damped harmonic oscillator (in the scenario i am looking for)? The idea of a simple damped harmonic oscillator was hard to understand at first, but driven seems to harder to learn, and to be honest, im at the point im stuck at because its something I want/need to learn for my game but there’s hardly any resources that explain it in a reasonable way.

Yes, I can’t seem to find resources on your specific application. On second thought, I’m not sure that the idea of a driven harmonic oscillator is useful in this context (my bad!). However, the problem remains the same; that is, accounting for the force induced by dragging one endpoint of the spring.

Actually, now that I’ve put this on pen and paper, the solution is much simpler. If you throw in the differential equation for a damped spring, you get what is essentially a PI controller:

a = w^2 * (X-x) - 2 * d * w * dx/dt
a is acceleration
w is undamped frequency
d is damping ratio
X is the target value

When you integrate this equation over time, you get the following:

v - v0 = w^2 ∫ (X-x) dt - 2*d*w * (x-x0)

or

v = B + w^2 ∫ E dt + 2*d*w * E
where
  E = X - x (the error with respect to the target)
  B = v0 - 2*d*w * (X-x0)

B is a bias term which accounts for the initial conditions. It is adjusted each time the setpoint is changed, to ensure that there aren’t any “bumps” in the velocity.

This is equivalent to a PI (Proportional, Integral) controller, which is why I suggest you learn about PID (Proportional, Integral, Derivative) controllers. They’re efficient, widely applicable, and in this case mathematically equivalent to what you’re trying to do. On the other hand, you could use the other velocity equation, but I feel obligated to inform you that something as useful as PID controllers exist.

Either way, you’ll want to take these equations and convert them into iterative algorithms. On that note, whenever you have some sort of differential equation involved, you can always discretize it to obtain an iterative solution. Computer programmers do this a lot.

7 Likes

Sorry for not being able to reply to this earlier, tbh I’ve been busy with school as I have exams coming up and I’ve lost alot of motivation within the last few days just because I can’t seem to get over this hurdle of knowledge I can’t seem to get/understand. However, I will try to sound as enthusiastic as possible :joy:

I’ve noticed how some fps games (atleast in roblox) have a spring like mechanic when you move the gun, this is what I would like to replicate. Here is an example of what I mean:

The springs act when the camera is moved or the player has a change in Y axis (e.g jumping)

Funny I literally gave up halfway through writing this im so out of it lately, tired af, stressed and out of motivation, hmmm :frowning:

Anyway, that’s the context to it, I’ll eventually see if I can get your code to work.

2 Likes

Sorry, just got around to actually looking at this properly, can you explain what each actual variables are and what they should be set to? You’ve explained some of them I just don’t want to play the guessing game of trying to figure out which one is which. I think these are the values you haven’t explained: dx, dt, v, v0, (idk what this symbol means either), x, x0, E, B.
Sorry if thats alot of values to ask description of but it would help me out alot.
Thanks again.

Let’s start off with the mathy stuff. ∫ is an integration symbol, and dx and dt are symbols that tell you which variable you’re integrating over (namely, x and t). These are concepts from calculus, so this may seem really scary. However, numerical integration is really easy to do with a computer program, so you needn’t worry. It’s basically just summing up the values of a function over a certain range. I would google “Numerical integration” and “Riemann integration”.

v is the velocity, and v0 is the initial velocity. Similarly, x is position and x0 is the initial position. E is the displacement of the target from the current position, i.e. x + E = X. It’s called ‘E’ because it’s the Error.

B is more subtle, it’s basically an extra term that you use to keep the velocity from jumping discontinuously when you change the target value. Every time you change the target, you have to adjust B to compensate for the sudden target change, so that the velocity doesn’t jump.

3 Likes

Sorry to ask about integrals, but I don’t learn about calculus till next year (A-levels) and tbh having to learn something just by reading wikipedia pages and other stuff is really hard to get a grip of understanding something.
I don’t mean to take time out of your day, but can you explain it to me? I’m basically completely new to this as we dont really learn much mathy stuff in my year (most we learnt was cos rule, circle theorems and sin rule), so stuff like this looks completely new I kinda need a teacher to explain it to me, so can you be my teacher for this?
Honestly I appreciate anything you tell me or teach me,I saw something explain it using functions like f(x) and g(x) but they introduce random variables such as D or n or something and I don’t understand where they come from or what they mean? If you could so kindly too put it into context for this?
I know I sound like i’m relying on you alot for this, but you have helped me alot so far :smiley:
Also quick question, how do you mean “initial” (e.g initial velocity, initial position)?
Thank you so much for your help.

I’ll start with the easy question first: the “initial” velocity is the value of the velocity at the instant we start simulation. This could be either when we start the program, or when we change the target value.

Integration, at least from a numerical standpoint, isn’t as hard as you think. Assuming E(t) is a function of time, we can calculate the integral of E like so:

local t, sum = 0
while t<1 do
    local dt = wait()
    sum = sum + E(t) * dt
    t = t + dt
end
print(sum)

Quite simple, practically speaking. Note that dt stands for delta t, or the change in t. If our variable were x, we would call it dx.

If you graph the curve of E(t), you can think of this as approximating the area under E(t) between 0 and 1, via a series of rectangles with height E(t) and length dt. This is called a Riemann sum, and this image from the Wikipedia page is instructive:

https://upload.wikimedia.org/wikipedia/commons/2/2a/Riemann_sum_convergence.png

In our example, the integration is from t=0 to t=1, but in the above image the integration starts at a negative value of t. It is important to know that, if E(t) is negative, the area is interpreted as negative.

3 Likes

I haven’t yet learnt any calculaus but I have tried to learn differentiation…
If we are differentiating are we doing ∫ (X-x) dt or ∫ (X-x) dt - 2*d*w * (x-x0)
I now know if you differentiate ax^b you get bax^b-1 but Im not too sure how to differentiate that…
Anyway not sure if I should be learning differentiation first or integration but I’m trying to understand…
Thank you again

For this specific application, you don’t need to know the derivative. Technically what you need is just a PI controller.

However, derivatives are also pretty easy to compute. You take the current value minus the last value, then divide by the time step. Example:

local t0 = tick()
local cur, last = E(0)
local dt = wait()
cur, last = E(tick()-t0), cur
print((cur-last)/dt)

And that will print the derivative (aka the rate of change) of E. There are more sophisticated ways of doing this which yield less noise-sensitive results, but for PID controllers I think this method is usually used.

3 Likes

I think I acctually managed to confuse myself here, unless you already said what is the actual function of ‘E’, or, how is it calculated, sorry for the slow reply i still really want to figure this out :smiley:

Persistence is a virtue for programmers! I’m glad you’re still working on it. I did mention the definition of E in an earlier post though:

Again, x is the current position, and X is the target position.

2 Likes

I think the reason I’m confused is when you said how to get the integral of X:

“Assuming E(t) is a function of time, we can calculate the integral of E like so:”

local t, sum = 0
while t<1 do
    local dt = wait()
    sum = sum + E(t) * dt
    t = t + dt
end
print(sum)

I’m not sure you actually said what E(t) actually was which is why I dont really understand.
Also why does t have to be less than 1? is there a reason?

Edit: Sorry for the slow reply again, lots of personal stuff going on and someone who i’m trying to fix a friendship with :frowning:

1 Like

I was slightly confused reading this, but I think I know what you’re confused about. In that case, let me clarify that I am using E and E(t) to mean the same thing. I guess I’ve been doing math for too long, so this might not make sense to a programmer. (When I was younger and more focused on programming, I had the opposite tendency: once in math class, I said sin(x) “returns” Opposite/Adjacent, and my math teacher insisted I say that sin(x) is Opposite/Adjacent. How the times change.)

I suppose that from a programmer’s perspective, E would either be a dynamic variable whose value at any given time is the error, or a function E() that returns the error when called. Does this clear things up?

No particular reason. In this case I integrated E between t=0 and t=1. For a PID controller, you’d update every frame, instead of integrating over a time window. Here’s an example:

local t, sum = 0
local dt
while true do
    dt = wait()
    sum = sum + E()*dt
    t = t+dt
    --
    print(sum)
end

You’d probably want to bind this to RenderStepped or some other event, so you have more control over the control flow, but you get the basic idea.

That’s tough. I admire your persistence in spite of that. Personally, I’m not a good multitasker, so something like that would take up my full attention. Good on you for keeping at it, and good luck with your friend.

3 Likes

Sorry, think im even more confused now :joy: How would you actually calculate the value of E? and if E is a function what would be the code inside of it? I’m not exactly sure what is meant by the error? Also, this might seem like a dumb question but what does sum = sum + E()*dt do? I know dt will be the time but im confused about why and what E()*dt is doing?

I guess it makes sense that its integrated however long it is, but would the total time passed be reset when the target moves?

Also thanks for the luck, im trying :slight_smile:

1 Like

Well, if E is such that x+E = X, then it follows that E = X-x, right? Then you should be able to calculate E from the target value and the current value.

Imagine a rectangle with height E and length dt. The area of this rectangle would be E*dt. It is negative when E is negative. Now plot E as a function of time. By adding a series of rectangles together, we can compute the (approximate) area between the curve of E and the time axis. Refer back to the picture I posted with the rectangles, if you need a visual.

2 Likes

Wow I have no idea how this went over my head! That sounds so easy, if I have any more questions I’ll definitely ask you about them.

That makes sense I guess, the visual definitely helped!
Yet again, thank you!

2 Likes

Okay, so far I think i’ve managed to create variables for everything, however I’m not a hundred percent sure they all of the right type.

local function create(c, A, B)
	local s = {}
	
	--a is the accelaration
	s.a = 0
	--w is the undamped frequency
	s.w = 0.5
	--d is the damping ratio
	s.d = 0.5
	--X is the target value
	s.X = Vector3.new()
	--x is the current value
	s.x = Vector3.new()
	--x0 is the initial position
	s.x0 = Vector3.new()
	--E is the error (calculated by doing X-x)
	s.E = 0
	--v0 is the initial velocity (set at the start or when target changes)
	s.v0 = Vector3.new()
	--v is the current velocity
	s.v = Vector3.new()
	--dx is the total time passed?
	s.dx = 0
	--dt is the time between each frame?
	s.dt = 0
	--B is the bias term (v0 - 2*d*w * (X-x0))
	s.B = 0
	
	
	s.update = function(x)
		
	end
	return s
end

Of course I still have the parameters c,A, and B where c is the position the spring is created and A and B are the variables that control the spring behaviour.

Also, I’m wondering why there is a complex solution to calculating v?
v = B + w^2 ∫ E dt + 2*d*w * E
What does it actually do or is it calculating?
Sorry for asking, I don’t like to code something without understanding it, because I made that mistake in the past and ended breaking everything and had to learn it anyway :joy: :man_shrugging:

Thank you again!! :smiley:

1 Like

a and dx don’t really belong, nor does dt since that should be local to the main loop. If dx were actually part of our setup, it would represent the change in position between the current frame and last frame. Nevertheless, we don’t need dx, so I digress.

Essentially, the second component 2*d*w * E responds immediately (it is Proportional) to the error, and it does most of the work. It is only dependent on the current error. The first term w^2 ∫ E dt can be thought of as summing up the past behavior, weighting it based on how long the error persists. Thus, it’s slower acting, but is useful for correcting errors that persist, especially steady-state error.

The variables d and w are derived from the differential equation for a damped harmonic oscillator. For a PID controller the equation would look different, but basically amounts to a change of variables.

At this point I feel like I’m basically translating Wikipedia pages into simple English, so I’d recommend some research. Some do seem to appreciate the information I’m giving, but a lot of this stuff has been explained on the internet better than I can. Here are some useful, non-technical resources I found:

https://www.csimn.com/CSI_pages/PIDforDummies.html

2 Likes