ODE45 is taking hours and hours to compute

I want to solve 27 odes, for that i formed equations with matrices.

First i formed A(27*27 matrix), B(27*27 matrix) and C(25*25 matrix). All three are interms of 'theta'(i used 'sym' for forming A, B, C matrices). Now iam going to use these matrices to form ode eq's and solving using ode45.

Here i gave script after forming A, B, C matrices for some confidentiality.

% A, B, C matrices formed interms of theta 

myfun = @(t,y)scriptname(t,y,A,B,C);

% dummy values for tspan and y0
tspan = [0 1];
y0 = zeros(27, 1);
% ode solver
sol = ode45(myfun,tspan,y0);

h = figure;
% plot 
plot(sol.x,sol.y(i,:)); 
    
function dydt = scriptname(t,y,A,B,C)
Wr = 2*pi*50;
p =2;
% evaluation of C (numerical) with theta = y(27)
Cn = double(subs(C,y(27)));
for i=1:25
    I(i,1)=y(i);
end
T1=1/2*p*I'*Cn*I

if t<0.5
    T2=0;
else
    T2=7.31;
end

V=[cos(Wr*t);
    cos(Wr*t+2.*pi/3.);
    cos(Wr*t-2.*pi/3.);
    zeros(21, 1);
    0;
    (T1-T2);
    y(26)]

% evaluation of A and B (numerical) with theta = y(27)
An = double(subs(A,y(27)));
Bn = double(subs(B,y(27)));
dydt = Bn\V-(An*y);
end

While running the script, it is taking hours and hours(i waited 5-6 hours and stopped compiling) but not giving any result.

I dont know what is wrong with the script.

Can anyone suggest me how to get result quickly.

Answer 1

[Walter Roberson]

Pay attention to the fact that the if statement was removed from the code, and that instead the run was split into two pieces that pass in different T2 values. The mathematics used for ode45() is such that if you use two different branches of an if statement in a single call to ode45(), then there is a good chance that your code is wrong, and that you need to stop the integration at the boundary and then resume integration from where you left off.

% A, B, C matrices formed in terms of theta 

commonvars = unique([symvar(A), symvar(B), symvar(C)]);  %probably just theta
Afun = matlabFunction(A, 'vars', commonvars);
Bfun = matlabFunction(B, 'vars', commonvars);
Cfun = matlabFunction(C, 'vars', commonvars);

tspan1 = [0, 0.5]; T2_1 = 0;
tspan2 = [0.5, 1]; T2_2 = 7.31;

myfun1 = @(t,y)scriptname(t, y, T2_1, Afun, Bfun, Cfun);
myfun2 = @(t,y)scriptname(t, y, T2_2, Afun, Bfun, Cfun);

y0_1 = zeros(27, 1);
% ode solver
[t_1, y1] = ode45(myfun1, tspan1, y0_1);
y0_2 = y1(end,:);
[t_2, y2] = ode45(myfun2, tspan2, y0_2);

t = [t_1; t_2];
y = [y1; y2];

h = figure;
% plot 
plot(t, y); 
    
function dydt = scriptname(t, y, T2, Afun, Bfun, Cfun)
    Wr = 2*pi*50;
    p =2;
    % evaluation of C (numerical) with theta = y(27)
    Cn = Cfun(y(27));
    for i=1:25
        I(i,1)=y(i);
    end
    T1=1/2*p*I'*Cn*I
    
    V=[cos(Wr*t);
        cos(Wr*t+2.*pi/3.);
        cos(Wr*t-2.*pi/3.);
        zeros(21, 1);
        0;
        (T1-T2);
        y(26)]
    
    % evaluation of A and B (numerical) with theta = y(27)
    An = Afun(y(27));
    Bn = Bfun(y(27));
    dydt = Bn\V-(An*y);
end

[Bathala Teja]

@Walter Roberson

thank you it worked, but waveforms are not as expected.

your steps for substituting y(27) are correct right??

[Walter Roberson]

My steps for using y(27) to calculate the A, B, C values are correct provided that your A, B, and C only have one symbolic variable in them.

[Bathala Teja]

I gave my entire script above

In that A = G, B = L, C = g.

Actually i used two variables theta and phi. But at the end of formation of G, L, g matrices "phi" will vanish out because of integration operation w.r.t phi.

So i will get G, L, g matrices interms of theta at the end.

And now i want to proceed to form and solve set of 1st order odes.

Usage of phi can be problem here???

[Walter Roberson]

I have made a number of small changes to my working copy.

The code has a number of integral() operations that cannot be expected to have closed-form solutions, but MATLAB has to try to integral each of them when it sees int()

You can substitute in vpaintegral() and get through to the creation of the G, L, g matrices.

However at that point you want to matlabFunction() . But matlabFunction() does not support vpaintegral() .

In a recent release, MATLAB added a 'hold' option to int() to get an unevaluated int() form that you could manipuate and later release() . However... matlabFunction() does not handle the 'hold' option either.

You end up having to upgrade the vpaintegral() objects into int() objects. Sadly, that is very time consuming.

You have to understand some of the more obscure parts of the Symbolic Toolbox in order to upgrade the vpaintgral() into int() objects. But doing so is mostly a waste of time as you can predict that very few of the expressions will turn out to have closed form solutions, but upgrading to int() requires that MATLAB spend a lot of time trying to do the integration.

It might perhaps be easier to find all of the vpaintegral() calls and generate matlabFunction for the integrand, emit a wrapper function that does integral(), and then replace the vpaintegral() with a symoblic function reference... I am not sure at the moment that all this can be made to work.

[Bathala Teja]

well i appreciate your explanation, but sorry to say that, i didnt understand.

I used symbolic math tool box for formulating my G, L, g matrices, is that a problem??

if it is the problem, i will try numercal way.

Do i need to replace int() with vpaintegral()??

Answer 2

[Jan]

Symbolic omputations need a lot of time. Can you implement the code numerically?

If the equation to be integrated is stiff, ODE45 tries to reduce the stepsize to mikroskopic values. Use a stiff solver in this case. e.g. ODE23S.

if t<0.5

Remember that Matlab's ODE integrators are designed to handle smooth functions only. Maybe this is a hard jump and the stepsize controller fails to pass this point. The correct way is to stop the integration at such jumps and restart it with the changed parameter.

[Bathala Teja]

i cant implement numerically since its difficult to form my A, B, C matrices numerically.

I also thought it will take somuch time but how much, i waited for 6-7 hours but still its running thats why i quit debugging.

I tried by removing these steps>>

if t<0.5

T1 = 0;

else

T1 = 7.31;

end

But no improvement same problem again.

Can this issue solved using ode23s???

[Walter Roberson]

i cant implement numerically since its difficult to form my A, B, C matrices numerically.

Outside of the ode call, use matlabFunction to generate function handles from A, B, C, and pass in the function handles and use those instead of subs()

[Bathala Teja]

means?

placing

An = Matlabfunction(A, theta, y(27));

Bn = Matlabfunction(B, theta, y(27));

Cn = Matlabfunction(C, theta, y(27));

lines before "sol = ode45......" line

??

[Stephen23]

@Bathala Teja: define those three (or however many) function handles before calling ODE45.

Then pass the function handles using the method shown here:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html#bsgprpq-8

[Bathala Teja]

sorry to bother you, i didnt understand your line "Outside of the ode call, use matlabFunction to generate function handles from A, B, C, and pass in the function handles and use those instead of subs()".

Can you show me in lines how can i replace with that double(subs())??

[Bathala Teja]

are you saying placing these below lines before ode45 line

An = matlabFunction(subs(A, y(27)))
Bn = matlabFunction(subs(B, y(27)))
Cn = matlabFunction(subs(C, y(27)))

It is not working...

[Steven Lord]

syms t
A = [1, t; t.^2, t.^3]

A = 

Afun = matlabFunction(A)

Afun = function_handle with value:
    @(t)reshape([1.0,t.^2,t,t.^3],[2,2])

B = subs(A, t, 2)

B = 

Bfun = Afun(2)

Bfun = 2×2

1 2 4 8

Generate Afun from A ahead of time outside your ODE function and pass it into your ODE function. When you would have substituted values into A instead evaluate Afun.

[Bathala Teja]

Have you understood my script which i mentioned above??

[Walter Roberson]

Yes, Steven and I have understood your script.

[Bathala Teja]

iam asking about my script>>Bathala Teja

If you understood can you edit my script to rectify the problem?

[Jan]

@Bathala Teja: The discussion would be much easier, if you post your complete code. Then the conversion to a direct numerical version is most likely easy. As long as the readers are guessing, they can give abstract hints only. If you do not understand these hints, a helpful answer is unlikely.

[Bathala Teja]

@Jan

any suggestions??