Physics prediction slightly off compared to physics simulation.

Physics
,
1

Hi all,

I am currently writing some physics prediction code (ballistic trajectory). I have time (flightDuration) and a distance (a simple float). I want to calculate the start velocity for a projectile.

I have the formula and it kinda works. Yet, the point of arrival is always a bit short and I don’t quite get why.

Used formula [1]:
9091582--1259299--upload_2023-6-20_10-10-10.png
Code: ```
var velocity = distanceVector / time - new Vector3(0, g, 0) * 0.5f * time;


I have created a single code example which you can add to an empty scene and run. It's really bare bones, yet it shows the problem. My expectation was that the white projectiles hit the red target exactly.

However, they are always a bit short. Like this (the white one should be at the position of the red target once it intersects the ground):

![9091582--1259311--upload_2023-6-20_10-20-34.png|967x546](upload://vz37qv7vMYERYz0DblMrV9Suac4.png)

[spoiler]

```csharp
using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class PhysicsPredictionTest : MonoBehaviour
{
    private Material targetMaterial;
    private Material groundMaterial;
    private Material projectileMaterial;
   
    IEnumerator Start()
    {
        var shader = Shader.Find("Standard");
       
        targetMaterial = new Material(shader);
        targetMaterial.color = Color.red;
       
        groundMaterial = new Material(shader);
        groundMaterial.color = new Color(0f, 0.7f, 0f);
       
        projectileMaterial = new Material(shader);
        projectileMaterial.color = Color.white;

        var ground = GameObject.CreatePrimitive(PrimitiveType.Quad);
        ground.name = "Ground";
        ground.GetComponent<MeshRenderer>().sharedMaterial = groundMaterial;
        ground.GetComponent<MeshCollider>().enabled = false;
        ground.transform.localRotation = Quaternion.Euler(90f, 0f, 0f);
        ground.transform.localScale = Vector3.one * 1000f;

        var target = GameObject.CreatePrimitive(PrimitiveType.Sphere);
        target.name = "Target";
        target.GetComponent<MeshRenderer>().sharedMaterial = targetMaterial;
        target.GetComponent<SphereCollider>().enabled = false;
        target.transform.localScale = Vector3.one * 0.05f;
       
        // Config (we derive the starting speed of the projectiles from these).
        float distance = 20f;
        float flightDuration = 3f;

        // Spawn target
        target.transform.localPosition = new Vector3(0f, 0f, distance);
       
        // Spawn projectiles
        for (int i = 0; i < 100; i++)
        {
            yield return new WaitForSeconds(0.1f);
            spawnProjectile(distance, flightDuration);
        }
    }

    void spawnProjectile(float distance, float time)
    {
        var startVelocity = calcStartVelocity(distance, time);
       
        var projectile = GameObject.CreatePrimitive(PrimitiveType.Sphere);
        projectile.name = "projectile";
        projectile.GetComponent<MeshRenderer>().sharedMaterial = projectileMaterial;
        projectile.transform.localScale = Vector3.one * 0.05f;
       
        var collider = projectile.GetComponent<SphereCollider>();
        collider.radius = 0.05f;
       
        var rb = projectile.AddComponent<Rigidbody>();
        rb.AddForce(startVelocity, ForceMode.VelocityChange);
    }
   
    Vector3 calcStartVelocity(float distance, float time)
    {
        var startPos = Vector3.zero;
        var targetPos = new Vector3(0f, 0f, distance);
        var distanceVector = targetPos - startPos;
       
        // With this it is always a bit short
        float g = Physics.gravity.y;
       
        // With this it is an exact match, even across vast distances (500+). But this does not compute with my brain.
        // float g = -9.87f;
       
        // Derived from equation for motion under constant acceleration.
        var velocity = distanceVector / time - new Vector3(0, g, 0) * 0.5f * time;
        return velocity;
    }
}

[/spoiler]

I expected my prediction to diverge from the simulation result on big distances due to floating point inaccuracy but right now it’s always off.

I feel like I am missing something very obvious (Unity 2021.2.0f1, 2022 LTS and 2023, no changes to physics settings)

9091582--1259314--upload_2023-6-20_10-23-11.png

Thank you.

Update: Just realized I posted it in the Scripting forum. If a mod is nearby maybe this should be moved to the physics forum. Thank you.

9091582–1259326–PhysicsPredictionTest.cs (3.08 KB)

2

I have a hunch now. If I put the starting position at f(t), where t = 0.5 * Time.fixedDeltaTime then I get reliable results.

Vector3 calcStartPosition(Vector3 startVelocity)
    {
        // Start at the position at half a physics update.
        float dT = Time.fixedDeltaTime * 0.5f;
     
        return new Vector3(
            0f,
            dT * startVelocity.y + Physics.gravity.y * dT * dT * 0.5f,
            dT * startVelocity.z
        );
    }

That got me thinking, maybe it’s a timing issue. After all, I was spawning a physics object in a coroutine and not FixedUpdate. Sadly, moving it to FixedUpdate did not change the outcome (it still requires the start-position-fix).

This is the current working version (spawning in FixedUpdate AND modifying the start position). If have tested it with distances up to 500 (working fine). Even tried different fixed delta times in the settings (up to 0.33). It also works fine.

Yet, I do not understand WHY it is working (or if it even should).

using System;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEngine.Serialization;

public class PhysicsPredictionTest : MonoBehaviour
{
    private Material targetMaterial;
    private Material groundMaterial;
    private Material projectileMaterial;
 
    public float Distance = 20f;
    public float FlightDuration = 3f;
 
    IEnumerator Start()
    {
        var shader = Shader.Find("Standard");
    
        targetMaterial = new Material(shader);
        targetMaterial.color = Color.red;
    
        groundMaterial = new Material(shader);
        groundMaterial.color = new Color(0f, 0.7f, 0f);
    
        projectileMaterial = new Material(shader);
        projectileMaterial.color = Color.white;

        var ground = GameObject.CreatePrimitive(PrimitiveType.Quad);
        ground.name = "Ground";
        ground.GetComponent<MeshRenderer>().sharedMaterial = groundMaterial;
        ground.GetComponent<MeshCollider>().enabled = false;
        ground.transform.localRotation = Quaternion.Euler(90f, 0f, 0f);
        ground.transform.localScale = Vector3.one * 1000f;

        var target = GameObject.CreatePrimitive(PrimitiveType.Sphere);
        target.name = "Target";
        target.GetComponent<MeshRenderer>().sharedMaterial = targetMaterial;
        target.GetComponent<SphereCollider>().enabled = false;
        target.transform.localScale = Vector3.one * 0.04f;
    
        // Spawn target
        target.transform.localPosition = new Vector3(0f, 0f, Distance);
    
        // Spawn projectiles
        for (int i = 0; i < 300; i++)
        {
            yield return new WaitForSeconds(0.1f);
            _spawnInNextUpdate = true;
        }
    }


    bool _spawnInNextUpdate = false;
 
    void FixedUpdate()
    {
        if (_spawnInNextUpdate)
        {
            _spawnInNextUpdate = false;
            spawnProjectile(Distance, FlightDuration);
        }
    }

    Vector3 calcStartPosition(Vector3 startVelocity)
    {
        // Start at the position at half a physics update.
        float dT = Time.fixedDeltaTime * 0.5f;
    
        return new Vector3(
            0f,
            dT * startVelocity.y + Physics.gravity.y * dT * dT * 0.5f,
            dT * startVelocity.z
        );
    }

    void spawnProjectile(float distance, float time)
    {
        var startVelocity = calcStartVelocity(distance, time);
    
        var projectile = GameObject.CreatePrimitive(PrimitiveType.Sphere);
        projectile.name = "projectile";
        projectile.GetComponent<MeshRenderer>().sharedMaterial = projectileMaterial;
        projectile.transform.localScale = Vector3.one * 0.05f;
        projectile.transform.localPosition = calcStartPosition(startVelocity); // New
    
        var collider = projectile.GetComponent<SphereCollider>();
        collider.radius = 0.05f;
    
        var rb = projectile.AddComponent<Rigidbody>();
        rb.AddForce(startVelocity, ForceMode.VelocityChange);
    }
 
    Vector3 calcStartVelocity(float distance, float time)
    {
        var startPos = Vector3.zero;
        var targetPos = new Vector3(0f, 0f, distance);
        var distanceVector = targetPos - startPos;
    
        // With this it is always a bit short
        float g = Physics.gravity.y;
    
        // Derived from equation for motion under constant acceleration.
        var velocity = distanceVector / time - new Vector3(0, g, 0) * 0.5f * time;
        return velocity;
    }
}
3

TL;DR; but I immediately spotted “time”. :wink:

I think you are on the right track. For one, using fixedDeltaTime is a must. But … how are you going to know how many fixedDeltaTime you need to calculate in advance?

Since fixedDeltaTime is a constant and probably 0.2 (I think) you will have some deviation for any time that is not a multiple of 0.2 like 1.19s may be a little short whereas 1.01s may noticably overshoot (or vice versa). I think you get what I’m thinking? Not sure if that is the issue but my gut says it is. :wink:

4

Thanks, yes, it may to have something to do with the physics simulation step. Though, I am pretty sure time or delta time should not have any effect on the outcome of the simulation. The shape of the parabola that is approximated should be the same in both scenarios (math and simlulation). At least to the point where it breaks apart due to float errors (long distance or big delta time).

The result of the physics simulation simply does not match the mathematical formula. The parabolas do not intersect the ground plane at the same position. The red target is where the mathematical formula hits the ground plane. The white projectile (green below) is where the simulation hits it.

In the screenshots in my first post you can see the white projectile nicely intersects the ground plane at y = 0. That, of course, does not happen with all distances as the physics step may hit a bit earlier or later. I chose the distance to get a nice intersection image. The thing is it should (if interpolated in between two steps) approximately hit the red target, yet it does not.

I expected the same outcode (same parabola with the math equation and the physics simulation). Yet I do get different ones. My current “solution” simply starts the parabola at another spot than the mathematical model and tbh, that bothers me as it’s not a logically sound fix but more a cover-up of the effects I see.

Here is a quick image of what I mean (I hope it explains it a bit better):
9092020--1259425--upload_2023-6-20_13-28-11.png

My current theory is that the physics engine already calculated one tick for the projectile when I add the force and therefore it is a bit displaced (downwards dueto gravity) causing the shortened parabola effect which I then counteract with the changing the start position (move upwards). Though I still have to confirm that (it’s just an idea at the moment). - I understand this may all be hard to follow if not familiar with the code. Thanks to anyone who follows me down into that rabbit hole :smile:

5

Turns out reducing the delta time fixes it too. I admit I have not yet tried that as it’s not my favourite solution (costly in terms of performance).

It seems the simulation simply diverges a lot even with the default delta time of 0.02f. I did not expect that. If set to 0.002f it all just works.

9092107--1259431--upload_2023-6-20_14-10-7.png

I guess my question has now become about how to adequatly predict the simulation divergence and counteract it so that it gives reasonably accurate results.

In my tests with dT = 0.02 if have seen it always shoots too short (never too far). Maybe that can be used to “fix” it. Hm, hmm …

Update:
I have now tried the “start-position-fix” with various Unity versions and distances and I can report (at least in the simple example posted above) it works very accurately up to at least a distance of 1000 units (hits objects with a size of 0.04). Whether or not this will hold true for other scenarios I can not say but it’s good enough for me. The neat thing is, since the fix is based on Time.fixedDeltaTime, it scales down automatically as delta time is reduced.

6

I’ve found a happy medium between performance and accuracy at quadruple an ascertainable target rate.

So if you want accuracy at 1/60th of a second “accuracy”, set fixedUpdate to 1/240th of a second.

7

Yes. Sadly that also means 4x the physics calculations are done, which may not be feasible in all scenarios.

8

Yeah… it’s a pain in the arse. There’s an OnAudioFilterRead that you can set to callback more frequently, but has complications of its own. However it’s in its own thread, so you can do as much as you like in there without causing too many issues to anything else.

Or, if you wanna go nuts, you can spin up your own thread and set a timer of your own.

9

For anyone interested. Here is the final “proof of concept” implementation. The code is a bit messy but it includes different solution strategies like calc start properties based on different known properties:

  • distance
  • speed
  • angle
  • duration
  • lowest used energy (this one I like very much as it gives you a nicely shooting turret without much thinking)

9093016--1259743--HittingTargetsWithBallisitics.gif
Code:

using System;
using System.Collections;
using UnityEngine;
using Quaternion = UnityEngine.Quaternion;
using Vector3 = UnityEngine.Vector3;

public class PhysicsPredictionTest : MonoBehaviour
{
    private Material targetMaterial;
    private Material groundMaterial;
    private Material projectileMaterial;
 
    public Transform Source;
    public Transform Target;
    public float FlightDuration = 2f;
    public float Angle = 45f;
    public float Speed = 20f;
 
    IEnumerator Start()
    {
        var shader = Shader.Find("Standard");
   
        targetMaterial = new Material(shader);
        targetMaterial.color = Color.red;
   
        groundMaterial = new Material(shader);
        groundMaterial.color = new Color(0f, 0.7f, 0f);
   
        projectileMaterial = new Material(shader);
        projectileMaterial.color = Color.white;

        var ground = GameObject.CreatePrimitive(PrimitiveType.Quad);
        ground.name = "Ground";
        ground.GetComponent<MeshRenderer>().sharedMaterial = groundMaterial;
        ground.GetComponent<MeshCollider>().enabled = false;
        ground.transform.localRotation = Quaternion.Euler(90f, 0f, 0f);
        ground.transform.localScale = Vector3.one * 1000f;

        var target = GameObject.CreatePrimitive(PrimitiveType.Sphere);
        target.name = "Target";
        target.GetComponent<MeshRenderer>().sharedMaterial = targetMaterial;
        target.GetComponent<SphereCollider>().enabled = false;
        target.transform.localScale = Vector3.one * 0.04f;
   
        // Spawn target
        target.transform.localPosition = Target.position;
   
        // Spawn projectiles
        for (int i = 0; i < 1900; i++)
        {
            yield return new WaitForSeconds(0.1f);
            _spawnInNextUpdate = true;
        }
    }

    IEnumerator deleteProjectile(GameObject go)
    {
        yield return new WaitForSeconds(FlightDuration + 1f);
        GameObject.Destroy(go);
    }


    bool _spawnInNextUpdate = false;
    Rigidbody _rb;
 
    void FixedUpdate()
    {
        if (_spawnInNextUpdate)
        {
            _spawnInNextUpdate = false;
            var rb = spawnProjectile(Source, Target, FlightDuration, Angle);
            if (_rb == null)
                _rb = rb;
            if(rb != null)
                StartCoroutine(deleteProjectile(rb.gameObject));
        }
    }

    protected bool useHighSolution = false;
 
    Rigidbody spawnProjectile(Transform source, Transform target, float time, float angle)
    {
        // By time
        //var startVelocity = calcStartVelocityByTime(source, target, time);

        // By angle
        /*
        bool possible = calcStartVelocityByAngle(out var startVelocity, source, target, angle);
        if (!possible)
            return null;
        */

        // By start speed
        /*
        useHighSolution = !useHighSolution;
        bool possible = calcStartVelocityByStartSpeed(out var startVelocity, source, target, Speed, useHighSolution);
        if (!possible)
            return null;
        */

        // By lowest energy needed
        var startVelocity = calcLowestStartVelocity(source, target);

        var projectile = GameObject.CreatePrimitive(PrimitiveType.Sphere);
        projectile.name = "projectile";
        projectile.GetComponent<MeshRenderer>().sharedMaterial = projectileMaterial;
        projectile.transform.localScale = Vector3.one * 0.1f;
        projectile.transform.localPosition = source.transform.position + calcStartPositionDelta(startVelocity); // New
   
        var collider = projectile.GetComponent<SphereCollider>();
        collider.radius = 0.1f;
   
        var rb = projectile.AddComponent<Rigidbody>();
        rb.AddForce(startVelocity, ForceMode.VelocityChange);

        return rb;
    }

    Vector3 calcStartPositionDelta(Vector3 startVelocity)
    {
        // Start at the position at half a physics update.
        float dT = Time.fixedDeltaTime * 0.5f;
   
        return new Vector3(
            0f,
            dT * startVelocity.y + Physics.gravity.y * dT * dT * 0.5f,
            dT * startVelocity.z
        );
    }
 
    Vector3 calcStartVelocityByTime(Transform source, Transform target, float time)
    {
        var path = target.position - source.position;
        var distanceXZ = new Vector3(path.x, 0f, path.z); // Project onto xz plane
        float distance = distanceXZ.magnitude;
        float altitude = path.y;
   
        // Solve in 2D (ZY plane)
        //   Derived from equation for motion under constant acceleration.
        //   var velocity = distanceVector / time - new Vector3(0, g, 0) * 0.5f * time;
        float g = Physics.gravity.y;
        var velocity = new Vector3(
            0f
            , altitude / time - g * time * 0.5f
            , distance / time
        );
   
        // Rotate back to 3D and return
        float angleY = Vector3.SignedAngle(distanceXZ, Vector3.forward, Vector3.up);
        return Quaternion.Euler(0f, -angleY, 0f) * velocity;
    }
 
    public bool calcStartVelocityByAngle(out Vector3 velocity, Transform source, Transform target, float angle)
    {
        if (source == null || target == null || float.IsNaN(angle))
            throw new Exception("Missing data source, target or angle.");

        var path = target.position - source.position;
        var distanceXZ = new Vector3(path.x, 0f, path.z);  // Project onto xz plane
        float distance = distanceXZ.magnitude;
        float altitude = path.y;
   
        // Solve in 2D (ZY plane)
        float g = Physics.gravity.y;
        var cos = Mathf.Cos(angle * Mathf.Deg2Rad);
        var value = distance * distance * -g /
                    (distance * Mathf.Sin(2f * angle * Mathf.Deg2Rad) - 2 * altitude * cos * cos);
        // If value < 0f then the target can not be reached with the given angle.
        if (value < 0)
        {
            velocity = Vector3.zero;
            return false;
        }
   
        var speed = Mathf.Sqrt(value);

        // Rotate back to 3D and return
        float angleY = Vector3.SignedAngle(distanceXZ, Vector3.forward, Vector3.up);
        velocity = Quaternion.Euler(-angle, -angleY, 0f) * new Vector3(0f, 0f, speed);

        return true;
    }
 
    /// <summary>
    /// This has two possible solutions. One with a high altitude (useHighSolution = true) and one with a low altitude (useHighSolution = false).
    /// </summary>
    /// <param name="startAngle"></param>
    /// <param name="source"></param>
    /// <param name="target"></param>
    /// <param name="speed"></param>
    /// <param name="useHighSolution"></param>
    /// <returns></returns>
    /// <exception cref="Exception"></exception>
    public bool calcStartVelocityByStartSpeed(out Vector3 startAngle, Transform source, Transform target, float speed, bool useHighSolution = false)
    {
        if (source == null || target == null || float.IsNaN(speed))
            throw new Exception("Missing data source, target or speed.");

        var path = target.position - source.position;
        var distanceXZ = new Vector3(path.x, 0f, path.z);  // Project onto xz plane
        float distance = distanceXZ.magnitude;
        float altitude = target.position.y - source.position.y;
   
        // Solve in 2D (ZY plane)
        var sqrSpeed = Mathf.Pow(speed, 2);
        var qadSpeed = Mathf.Pow(speed, 4);
        var x = distance;
        var y = path.y;
        float g = Physics.gravity.y;
        var sqrValue = qadSpeed - g * (g * x * x + 2 * -y * sqrSpeed); // Not sure why -y (wikipedia says +y but hey, it works).
   
        // If value < 0f then the target can not be reached with the given angle.
        if (sqrValue < 0f)
        {
            startAngle = Vector3.zero;
            return false;
        }
   
        // Has two solutions, see: https://en.wikipedia.org/wiki/Projectile_motion
        var angleA = Mathf.Atan((sqrSpeed + Mathf.Sqrt(sqrValue)) / (g * x)) * Mathf.Rad2Deg;
        var angleB = Mathf.Atan((sqrSpeed - Mathf.Sqrt(sqrValue)) / (g * x)) * Mathf.Rad2Deg;
   
        // Rotate back to 3D and return
        float angleY = Vector3.SignedAngle(distanceXZ, Vector3.forward, Vector3.up);
        startAngle = Quaternion.Euler(useHighSolution ? angleA : angleB, -angleY, 0f) * new Vector3(0f, 0f, speed);
        return true;
    }
 
    /// <summary>
    /// Find the angle with the lowest needed start speed and use that.
    /// </summary>
    /// <param name="velocity"></param>
    /// <param name="source"></param>
    /// <param name="target"></param>
    /// <returns></returns>
    Vector3 calcLowestStartVelocity(Transform source, Transform target)
    {
        // Wikipedia: The launch should be at the angle halfway between the target and Zenith (vector opposite to Gravity)
        // See: https://en.wikipedia.org/wiki/Projectile_motion
        var path = target.position - source.position;
        var angle = Vector3.Angle(path, Physics.gravity) / 2f;
        calcStartVelocityByAngle(out var velocity, source, target, angle);
        return velocity;
    }
}

May it be useful to all who come after :slight_smile:

10

Game physiscs engines don’t use the exact formulas, they incrementally update the solution with some timestep.

I am slightly simplifying thing but roughly speaking every physics engine step does something like this.

f(n) = gravity + other forces
a(n) =  f(n)/m
v(n)=v(n-1)  + a(n)* dt
position(n) = position(n-1)+ v(n) * dt

Assuming you had any calculus class, even if you were half a sleep, you probably remember drawings of approximatimg integral by drawing a bunch of rectangles. Physics engine is doing essentially the same, every update it draws a new rectangle to approximate the integral of motion equation. But that’s that - an approximation. It doesn’t produce the same result as exact formula even if you ignore floating point stuff. If you reduce the timestep (draw more thinner rectangles) result you get a better approximation which is closer to the exact formula, but in general case you don’t get exact result.

There is even whole sub-field/problem of math studying how to better approximate functions by incrementally calculating them. Read Numerical methods for ordinary differential equations - Wikipedia for more details. Different techniques produce slightly different results, and some of them can produce more accurate approximations with less steps (potentially at the cost of each step being more complex). For example when calculating position(n) do you use v(n), v(n-1), (v(n)+v(n-1))/2 or something more complex.

Not every math problem has formula for exact solution, so sometimes such incremental approaches are necessary to do any kind of calculations. Also for real world engineerring problems, objects often have complex shapes and many different forces interacting, making it difficult to apply the exact formulas directly. In such situation only solution is splitting things into smaller chunks and timesteps until you get a good enough approximation.

Situation with game physics engines is slightly different. Simple motion equations have the formulas for exact solutions, but those are only usable in simple textbook situations where you know ahead of time everything that will happen. Games are typically interactive, so at any point of time player can potentially perform an action changing the conditions of simulation, other problem that in complex environments with many objects bouncing around you don’t know what will collide with what and when. Unlike simple ball thrown in empty space where only that can change it’s trajectory is hitting ground.

A phyisics engine needs to choose a tradeoff (and corresponding methods) whether it’s faster to do more simpler steps, or less but more complicated steps. The optimal point on real hardware may vary depending on microarchitectural details of specific cpu.

Number of steps and compuation speed isn’t only tradeoff, there is also accuracy and stability.

For example if you simulate an object attached to a spring, it should keep oscilating with fixed amplitude, and if there is little bit of damping it after oscillating for a while should settle. Even if the behavior doesn’t perfectly match with exact formula you would probably want your simulation to have similar behavior. But if the numerical method is chosen or implemented badly you may endup in situation, where with large simulation step instead of settling, integration error accumulates, and bigger the error faster the error keeps growing until physics simulation explodes.
Here is a website https://myphysicslab.com/ myPhysicsLab Single Spring for simple physics simulations which allows you to experiment with different timesteps and numerical methods . With the default settings you should see the spring settle. But if you increase the timestep and switch to Eulers method, you might observe that instead of settling the amplitude increases.

For a specific situation if you fully expand and add up the formulas for an incremental approach, you will likely result with a formula that’s very similar to the exact solution but has small additional term on the scale of single timestep. If you know which exact numerical method a physics engine is using, you can cancel out this factor which is what you are trying to do above with your offset (or you can randomly guess until it matches). But that only works under very specific conditions, and if the situation slightly changes (for example there is ceiling and ball bounces once, wind resistance get introduced) not only you would have to change the formula for exact solution, you would also likely need a different formula for the adjustment factor.

It might not even be considered an API breaking change if physics engine update changed the numeric method it internally uses, but it would result in requiring a different adjustment factor to compensate between exact solution and the physics engine approximation. Although in practice a physics engine probably wouldn’t change it carelessly since people tend to depend on behavior not changing even if it’s result of nondocumented internal design decisions.