How to set up solving of multiple ODEs in Matlab propperly?

Question

I have a task of writing an M script that would set up (and finally solve) a set of 19 differential equations and all of the equation coefficients. I am not sure what is the best way to input those equations.

Example of some equations:

Example of coefficients:

I don't think that this would be suited for Simulink. Other examples use a function in the form of @(t,x) where t is time and x the vector of all variables.

There are loads of "examples" online, but they don't seem to be appropriate for such a large set of large equations.

For example this exmple of solving 3 equatons

Even for 3 simple equations as they have solved, the functions are getting messy:

f = @(t,x) [-x(1)+3*x(3);-x(2)+2*x(3);x(1)^2-2*x(3)];

Using this notation and getting to the x(19) and cross-referencing all instances of x would be a mess.

I would like your help, and a simple example of how could I write these equations line by line, maybe using the symbolic toolbox, and then put them in an array that I can then forward to the solver.

As I said, I know there are examples online, but they tackle only the most basic system, and really they don't scale well if you want a clean and easily readable code.

I would like to have a similar to Wolfram alpha, where you type variable names as they are (instead od x(1), x(2), ... m x(19)) and if it's possible to get all solution vectors with their variable names.

score 2 · Answer 1 · answered Sep 08 '19 at 19:01

You don't have to use an anonymous function handle as an ode function, you can create a separate function file (as shown in the odefun section of ode45.

For example, your odefun can look like:

function dy = myode(t,y)
    % first unpack state variables
    i_d2 = y(1);
    i_q2 = y(2);
    ...
    gamma2 = y(end-1);
    omega2 = y(end);

    % determine all constants
    c34 = expression_for_c34;
    ...
    c61 = expression_for_61;

    % determine state derivative
    i_d2_dot = expression;
    ...
    omega2_dot = expression;

    % pack state derivative based on order in state vector
    dy(1) = i_d2_dot;
    ...
    dy(end) = omega2_dot;

end

From this myode function you can also call other functions to e.g. determine the value for some coefficients based on the current state. Next, integrate the system using a suitable ode solver:

[t,y] = ode45(@myode,tspan,y0);

How to set up solving of multiple ODEs in Matlab propperly?

1 Answers1