Spring animation using Runge Kutta ?

Started by
3 comments, last by h4tt3n 15 years, 10 months ago
I'm trying to create a simplified animation of the one found here: http://www.myphysicslab.com/spring2d.html the only difference is my anchor is static at (200,200) and the ball is started at 50,50 and cannot be moved by the user. my animation is jacked...can someone help??.... please.. i currently have this hacked to look somewhat correct but i would like the correct solution using runge kutta also when should i use the runge kutta method? since i know i can create this effect without runge kutta... i could just setup a vector class and scale the vector accordingly manually without having to go into complicated differential equations. It seems like the only reason its used in that animation is to help in generating the graph - the transient response graph of any system. which would be difficult to do if i went the vector route. Is this the distinguishing factor on when to use runge kutta?

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Data;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Imaging;
using System.Windows.Navigation;
using System.Windows.Shapes;


namespace FirstAnimation
{
    public partial class Window1 : Window
    {   
        double time = 0;
        double dt = 0.03;
            // dt is constant to create a frame based dependence
            // time based dependence would need: dt = (renderArgs.RenderingTime – lastRender).TotalSeconds;


        const double HEAD_MASS = .15;
        const double GRAVITY_CONSTANT = 9.8;
        const double SPRING_REST_LENGTH = 25;
        const double SPRING_CONSTANT_K = 6;
        const double SPRING_DAMPING_CONSTANT = 0.25;


        // modified after eq solved @ dt
        private double Angle = 45 * Math.PI / 180;
        private double SpringStretch = 15;


        // dependent vars
        private double dependentPositionX = 50;
        private double dependentPositionY = 50;
        private double dependentVelocityX = 0;
        private double dependentVelocityY = 0;
        private double[] dependentVariables;



        public Window1()
        {
            InitializeComponent();

            this.time = 0;


            this.dependentVariables = new double[]
            { 
                dependentPositionX,
                dependentVelocityX,
                dependentPositionY,
                dependentVelocityY
            };
            CompositionTarget.Rendering += StartAnimation;
        }


        private void StartAnimation(object sender, EventArgs e)
        {
            ODESolver.Function[] fn = new ODESolver.Function[] 
            { 
                function_DUx, function_DVx,
                function_DUy, function_DVy
            };
            double[] result = ODESolver.RungeKutta4(fn, this.dependentVariables, this.time, this.dt);
            this.dependentVariables = result;           
            // this.time += this.dt;



            Point anchor = new Point(200, 200);
            Point currentBobblePosition = new Point(result[0], result[1]);
            double SpringLength_L = Math.Sqrt( Math.Pow((currentBobblePosition.X-anchor.X), 2) + Math.Pow((currentBobblePosition.Y-anchor.Y), 2));
                // Pythagorean theorem to find SpringLength_L


            this.SpringStretch = SpringLength_L - Window1.SPRING_REST_LENGTH;
            this.Angle = Math.Sinh((currentBobblePosition.X - anchor.X) / SpringLength_L);
            

            this.rectangle1.SetValue(Canvas.LeftProperty, currentBobblePosition.X );
            this.rectangle1.SetValue(Canvas.TopProperty, this.SpringStretch * Math.Cos(this.Angle));

            Console.Write("x = {0};  y = {1};  s = {2};  2 - {3}  3 - {4} \n", result[0], result[1], this.SpringStretch, result[2], result[3]);
        }






     
        //θ = angle (0 = vertical, increases counter-clockwise)
        //S = spring stretch (displacement from rest length)
        //L = length of spring
        //u = position of bob
        //v = u'= velocity of bob
        //a = u''= acceleration of bob
        

        private double function_DUx(double[] x, double t)
        {
            return x[1]; // Vx
        }
        private double function_DUy(double[] x, double t)
        {
            return x[0]; // Vy
        }
        private double function_DVx(double[] x, double t)
        {
            //x[0] == Vy
            //x[1] == Vx
            return (-1 * (Window1.SPRING_CONSTANT_K / Window1.HEAD_MASS) * this.SpringStretch * Math.Sin(this.Angle) - (Window1.SPRING_DAMPING_CONSTANT / Window1.HEAD_MASS) * x[1]);
        }
        private double function_DVy(double[] x, double t)
        {
            return (Window1.GRAVITY_CONSTANT - (Window1.SPRING_CONSTANT_K / Window1.HEAD_MASS) * this.SpringStretch * Math.Cos(this.Angle) - (Window1.SPRING_DAMPING_CONSTANT / Window1.HEAD_MASS) * x[0]);
        }
    }
}

following class was taken from Practical WPF Graphics Programming by Jack Xu


using System;
using System.Windows;
namespace FirstAnimation
{
    public class ODESolver
    {
        public delegate double Function(double[] x, double t);
        public static double[] RungeKutta4(Function[] f, double[] x0, double t0, double dt)
        {
            int n = x0.Length;
            double[] k1 = new double[n];
            double[] k2 = new double[n];
            double[] k3 = new double[n];
            double[] k4 = new double[n];

            double t = t0;
            double[] x1 = new double[n];
            double[] x = x0;

            for (int i = 0; i < n; i++)
                k1 = dt * f(x, t);

            for (int i = 0; i < n; i++)
                x1 = x + k1 / 2;

            for (int i = 0; i < n; i++)
                k2 = dt * f(x1, t + dt / 2);

            for (int i = 0; i < n; i++)
                x1 = x + k2 / 2;

            for (int i = 0; i < n; i++)
                k3 = dt * f(x1, t + dt / 2);

            for (int i = 0; i < n; i++)
                x1 = x + k3;

            for (int i = 0; i < n; i++)
                k4 = dt * f(x1, t + dt);

            for (int i = 0; i < n; i++)
                x += (k1 + 2 * k2 + 2 * k3 + k4) / 6;

            return x;
        }
    }
}

Advertisement
I haven't looked at the formulas in detail, but Math.Sinh() seems quite out of place... did you mean Math.Asin()?
I'd use Math.Atan2(), though; the pendulum might swing above horizontal.

Omae Wa Mou Shindeiru

Yup, sinh is incorrect, asin would be more appropriate. Atan2 not so much since he appears to have an edge length and a hypoteneuse, not two edge lengths. But, I don't have time to evaluate such a long piece of code for overall correctness, I'm afraid.

Am I correct in assuming this is C# code?
Graham Rhodes Moderator, Math & Physics forum @ gamedev.net
yes the language is c# code. its a wpf application.

fair enough, i knew the code snippet wouldn't be looked at much. Can you all answer the second part of my question though:

Quote:

also when should i use the runge kutta method?
since i know i can create this effect without runge kutta... i could just setup a vector class and scale the vector accordingly manually without having to go into complicated differential equations.


It seems like the only reason its used in that animation is to help in generating the graph - the transient response graph of any system. which would be difficult to do if i went the vector route. Is this the distinguishing factor on when to use runge kutta?



Well, The 4th order Runge-Kutta is simply a very good integration algorithm.

Anyway, if you want to make a spring simulation I recommend that you solve it numerically using force-, acceleration-, velocity-, and position vectors. First calculate the spring force on the bob using Hookes law of elasticity, then find its acceleration. Finally use RK4 to integrate with time and find velocity and position. This way you can simulate any combination of springs and point masses with the same - fairly simple - piece of generic code. Simulations based on trig functions doesn't have the same degree of freedom as vector based ones.

Cheers, Michael

This topic is closed to new replies.

Advertisement