Using ODEsolver to solve drag equation


#1

Hi!

I’m very new to C# and Math.Net programming. I’ve got some experience in other languages, but have mainly used MATLAB in recent years.

I’m trying to create a method that calculates the trajectory of a projectile under the influence of drag and gravity.
In addition the drag coefficient is velocity dependent (simple table lookup).

Let x be the horizontal distance and y be the vertical distance. The initial position, velocity and direction is known.
The forces acting on the projectile is:
F = m*g - rho * v^2 *A * CD/2
where rho is the density of the air, A is the presented area of the projectile and CD is the drag coefficient, all known.

How would I use the ODEsolvers in Math.Net to find y = f(x)?
The documentation is really lacking when it comes to the ODEsolvers…

I assume I need to create a method/function for the equation of motion and pass it to the solver?
How do I do that? Remember I’m really new to C#…

Hope you can help me!


(Peter Vanderwaart) #2

Differential equations are not my strong suit, but let me give it a try. It looks to me like the ODE solvers RungeKutta and AdamsBashforth are meant for functions of a single variable (or state). Your problem has four variables: x-position, y-position, x-velocity, y-velocity. They are not going to do what you need them to do.

Note: you may want to express your velocity vector in polar coordinates rather than (x,y) terms.


(Peter Vanderwaart) #3

Since I’d never used an ODE solver, and since the OP asked for some instruction, I thought I’d experiment. I thought it would be interesting to see how accurately the ODEsolvers would estimate a sine function over a number of cycles.

Program here:

    // use ODE solvers to estimate sine function
    // 10 cycles, i.e t from 0 to 20 pi

    // starting conditions
    double y0 = 0.0;
    double start = 0.0;  
    double end = 20.0 * Math.PI;
    int N = 150;
    var f1 = new Func<double, double, double>((t, y) => Math.Cos(t));

    var r = RungeKutta.SecondOrder(y0, start, end, N, f1);
    var s = AdamsBashforth.SecondOrder(y0, start, end, N, f1);

    MyChart.Visible = true;
    MyChart.Palette = Charting.ChartColorPalette.Bright;

    Charting.Title chtTitle = new System.Windows.Forms.DataVisualization.Charting.Title();
    Font chtFont = new System.Drawing.Font("Arial", 16);
    chtTitle.Font = chtFont;
    chtTitle.Text = "ODE Solver Examples";
    MyChart.Titles.Add(chtTitle);

    // Create Series 
    MyChart.Series.Clear();
    MyChart.Series.Add("Runge Kutta");
    MyChart.Series["Runge Kutta"].ChartType = Charting.SeriesChartType.Line; ;

    // populate Series    
    for (int i = 0; i < r.Length; i++)
    {
        double x = Convert.ToDouble(i) * 20.0 * Math.PI / Convert.ToDouble(s.Length);
        MyChart.Series["Runge Kutta"].Points.AddXY(x, r[i]);
    }

    MyChart.Series.Add("Adams-Bashforth");
    MyChart.Series["Adams-Bashforth"].ChartType = Charting.SeriesChartType.Line; ;

    // populate Series    
    for (int i = 0; i < s.Length; i++)
    {
        double x = Convert.ToDouble(i) * 20.0 * Math.PI/Convert.ToDouble(s.Length);
        MyChart.Series["Adams-Bashforth"].Points.AddXY(x, s[i]);
    }
    
 MyChart.SaveImage(@"C:\users\peter\documents\ODESolvers.jpg",Charting.ChartImageFormat.Jpeg);

The graph is here. You can see that using 15 points per cycle (about 24 degrees per point), the AB solver overestimates the upper and lower bounds and the KE underestimates them. If you triple the number of observations, the results seem perfect to the naked eye. I’m impressed that the results stayed so true over 10 complete cycles. It can be argued, and it’s probably true, that the sine function is not a difficult challenge compared to some other things.