How to achieve symbolic numerical integration

Hi Community,

I am trying to optimise a simulation so that it does not take 20 seconds to run a 3 second simulation. I will explain the simplified problem first then explain its context at the end (the function I am numerically integrating is far more complex and Matlab has failed to find an analytical solution to it).

Say I have a simple function, that must be integrated w.r.t. chi each iteration of ode45, such as:

funct_1 = @(chi, x) chi*x

The issue with doing a numerical integration each iteration of ode45, like so (where x is redefined each iteration od the ode solver and is the state variable):

integral(@(chi) fun(chi,x), 0, 1)

is that the majority of the simulation time is spent numerically integrating this function.

I have noticed that it is unneccessary to numerically intergate this each iteration as the bounds of chi never change, and therefore the numerical integration is really returning a function of x, where x is redefined each iteration.

So the hypothesis is that i can improve computation time by defining the numerical intergation of funct_1 w.r.t chi from 0 to 1 as a function itself of the state variable. Therefore the numeric integration is computed once and then values are simply plugged in.

i.e. I want to perform a numeric integration of funct_1 from chi = 0 to chi = 1, while keeping x symbolic (or a function variable).

How do I achieve this?

I have tried some solutions, but they haven't achieved my aims, such as the following:

funct_2 = @(x) integral(@(chi) fun(chi,x), 0, 1)

This returns:

funct_2 =

  function_handle with value:

    @(x)integral(@(chi)fun(chi,x),0,1)

Which is a function of x, but is really just a function of a function and is calling integral each iteration, and so has no impact on computation speed (in fact making it marginally slower).

I do not want the function to call a function, but instead for the function to be a function of x within it's own right (and to still exist if I were to delete the variable funct_1 for instance).

For context I am solving a physical dynamical system, where I have a generalised force of drag which is a function of the state variables and also the position on the body which is being acted on by air resistance. This generalised force of drag must be integrated across the whole body to get the generalised force of drag acting on the system, however the equation is too complex to be solved analytically and so must be solved numerically.

Numerically solving for the generalised force of drag each iteration takes very long, and in fact takes ~70x as long as the case that I would just assume a single force of drag acting at the centre of mass (and thereby foregoing the need to numerically intergate).

Any help would be appreciated.

[Dyuman Joshi]

" I want to perform a numeric integration of funct_1 from chi = 0 to chi = 1, while keeping x symbolic (or a function variable)."

You can not. You can either keep the integration numeric or symbolic, not a mix of both.

Though you can do symbolic integration w.r.t a single variable -

syms x y
F = x^2*y;
int(F, x, 0, 1)

ans = 

However, depending upon the complexity of the equation, that might not be as swift as you want to obtain.

What components does the drag force involve? Both - Pressure (including air resistance) and Viscous drag? Or is air resistance factored separately?

Could you provide additional details regarding the system?

[Matthew Ediz Beadman]

It is a simple system with the bare bones drag equation of

F_drag = -1/2 sqrt( dot(v_x,v_x) ) * rho * Cd * A * abs( cos(phi) ) * v_x

This is just the linear drag equation where -sign(v)*v^2 is switched for -norm(v)*dot(v,v) to make the equation a vector, and area is A*abs( cos(phi) ) due to how the system is formulated.

It uses Lagrangian Mechanics, so this is converted to a generalised force, and v_x is the velocity of a generic position 'x' on the rod. The idea is then to integrate the generalised force resulting from drag at the point x on the rod w.r.t x ranging from 0 (where the rod start) to l (which is a defined numeric parameter and signifies the rods end).

The generalised force of the rod is not a simple equation, and does not respond to analytical techniques of integration.

% drag coefficient
drag_coefficient = (-1/2).*Aa.*Cd.*rho.*abs(cos(psi_));

% generalised force of drag acting at a generic point on the rod
inf_drag_sym = drag_coefficient.*...
    [(x_d+psi_d.*chi.*cos(psi_)+(-1).*h2.*theta_d.*cos(theta)).*((x_d+psi_d.*chi.*cos( ...
      psi_)+(-1).*h2.*theta_d.*cos(theta)).^2+((-1).*psi_d.*chi.*sin(psi_)+h2.* ...
      theta_d.*sin(theta)).^2).^(1/2);0;(-1).*h2.*cos(theta).*(x_d+psi_d.*chi.*cos( ...
      psi_)+(-1).*h2.*theta_d.*cos(theta)).*((x_d+psi_d.*chi.*cos(psi_)+(-1).*h2.* ...
      theta_d.*cos(theta)).^2+((-1).*psi_d.*chi.*sin(psi_)+h2.*theta_d.*sin(theta)) ...
      .^2).^(1/2);chi.*cos(psi_).*(x_d+psi_d.*chi.*cos(psi_)+(-1).*h2.*theta_d.*cos( ...
      theta)).*((x_d+psi_d.*chi.*cos(psi_)+(-1).*h2.*theta_d.*cos(theta)).^2+((-1).* ...
      psi_d.*chi.*sin(psi_)+h2.*theta_d.*sin(theta)).^2).^(1/2)];

Answer 1

[Star Strider]

Without knowing more about ‘chi’ and the characteristics of its integrated result, one option (especially if the integration limits are always the same) could be to calculate it as a vector or matrix and then use the appropriate interpolation function (that are generally fast) to return the desired value. The approach would be similar to that described in ODE with Time-Dependent Terms.

If ‘chi’ and its integral are well-behaved (without discontinuities or singularities), this could work.

[Matthew Ediz Beadman]

That did the trick! Brought down the simulation time about 70%. I would have never thought of data interpolation, thank you :)

[Star Strider]

As always, my pleasure!