Out of memory - dsolve 3rd order system

I'm trying to use dsolve to solve the 3rd order LTI system with input and all initial conditions 0. When I run my code, it seems to run forever until I get an out of memory error (on desktop MATLAB 2021a) or my session expires (MATLAB Online). I can't figure out what the problem is. I based this code on similar examples that run fine, which were found here. I would greatly appreciate any help.

clear;
syms y(t) t;
Dy = diff(y, 1);
D2y = diff(y, 2);
f = 10*exp(-t);
sysB = diff(y,3) + 7*diff(y,2) + 8*diff(y,1) + 9*y == diff(f,3) + 2*diff(f,2) + 3*diff(f, 1) + 4*f;
cond1 = y(0) == 0;
cond2 = Dy(0) == 0;
cond3 = D2y(0) == 0;
conds = [cond1; cond2; cond3];
y1 = dsolve(sysB, conds)
ezplot(y1, 0, 10)

Errors:

Error using mupadengine/feval_internal

Out of memory.

Error in dsolve>mupadDsolve (line 334)

T = feval_internal(symengine,'symobj::dsolve',sys,x,options);

Error in dsolve (line 203)

sol = mupadDsolve(args, options);

Error in exam2_3 (line 14)

y1 = dsolve(sysB, conds)

render process terminated: null

java.lang.NoSuchFieldError: TS_PROCESS_OOM

Answer 1

[Torsten]

Does this one work better ?

syms t y1(t) y2(t) y3(t)
f = 10*exp(-t);
eqn1 = diff(y1) == y2;
eqn2 = diff(y2) == y3;
eqn3 = diff(y3) == -7*y3 - 8*y2 - 9*y1 + diff(f,3) + 2*diff(f,2) + 3*diff(f, 1) + 4*f;
eqns = [eqn1,eqn2,eqn3];
cond1 = y1(0) == 0;
cond2 = y2(0) == 0;
cond3 = y3(0) == 0;
conds = [cond1,cond2,cond3];
y = dsolve(eqns,conds,'MaxDegree',3)

[Chris G.]

I ran this for 2.5 hours without getting a result. It did not run out of memory though. Is a runtime longer than 2.5 hours expected?

[Torsten]

I think using pencil and paper would be faster ...

Or a numerical solution.

Answer 2

[Walter Roberson]

syms y(t) t;
Dy = diff(y, 1);
D2y = diff(y, 2);
f = 10*exp(-t);
sysB = diff(y,3) + 7*diff(y,2) + 8*diff(y,1) + 9*y == diff(f,3) + 2*diff(f,2) + 3*diff(f, 1) + 4*f;
cond1 = y(0) == 0;
cond2 = Dy(0) == 0;
cond3 = D2y(0) == 0;
conds = [cond1; cond2; cond3];
y1general = dsolve(sysB);

Dy1general = diff(y1general, t);
D2y1general = diff(Dy1general, t);

cond1general = subs(y1general, t, 0);
cond2general = subs(Dy1general, t, 0);
cond3general = subs(D2y1general, t, 0);

%the three cond*general are linear equations in the constants of
%integrations, so in theory they are trivial to solve. 
%in practice because the equations are long, it takes a fair time
%to calculate the coefficients. We can help some by guiding the order
%of solution. But even so this is slower than the dsolve() itself.

condeqn1 = [cond1general, cond2general, cond3general];
vars = symvar(condeqn1)

var1partial = solve(condeqn1(1), vars(1));
condeqn2 = subs(condeqn1(2:end), vars(1), var1partial);
var2partial = solve(condeqn2(1), vars(2));
condeqn3 = subs(condeqn2(2:end), vars(2), var2partial);
var3partial = solve(condeqn3, vars(3));

%and then you would do the back substitutions to recover the full constants

It might be faster to extract the coefficients of the linear conditions equations using coeff() and develop "model" equations like

p1 * C1 + q1 * C2 + constant1 == 0
p2 * C1 + q2 * C2 + r2 * C3 + constant2 == 0
p3 * C1 + q3 * C2 + r3 * C3 + constant3 == 0

and solve the model for C1, C2, C3 in terms of p1, p2, p3, q1, q2, q3, r2, r3... and then to subs() the actual coeffs in to the solution. 'Cuz the coefficients really are pretty long and doing all of the multiplications as you go takes a fairly long time... much faster to construct the general form and substitute values in.

By the way: the solutions involve complex values. For example the first condition turns out to be

C1 + C2 - 0.73490750273865887398095800213633 + 0.3174874795997589233399663703506i == 0

[Walter Roberson]

Note: the final steps building up the partials has been running for about 5 hours on my system.

Maybe I really should use the coeffs() approach, or at least equationsToMatrix

[Walter Roberson]

syms y(t) t;
Dy = diff(y, 1);
D2y = diff(y, 2);
f = 10*exp(-t);
sysB = diff(y,3) + 7*diff(y,2) + 8*diff(y,1) + 9*y == diff(f,3) + 2*diff(f,2) + 3*diff(f, 1) + 4*f;
cond1 = y(0) == 0;
cond2 = Dy(0) == 0;
cond3 = D2y(0) == 0;
conds = [cond1; cond2; cond3];
y1general = dsolve(sysB);

Dy1general = diff(y1general, t);
D2y1general = diff(Dy1general, t);

cond1general = subs(y1general, t, 0);
cond2general = subs(Dy1general, t, 0);
cond3general = subs(D2y1general, t, 0);

%the three cond*general are linear equations in the constants of
%integrations, so in theory they are trivial to solve. 
%in practice because the equations are long, it takes a fair time
%to calculate the coefficients. We can help some by guiding the order
%of solution. But even so this is slower than the dsolve() itself.

condeqn1 = [cond1general, cond2general, cond3general];
vars = symvar(condeqn1);

[A, b] = equationsToMatrix(condeqn1);
value_of_constants = double(A)\double(b)

value_of_constants = 
   0.0337 + 0.0713i
   0.7012 - 0.3888i
  -0.2694 - 1.3386i

y1specific = vpa( subs(y1general, vars(:), sym(value_of_constants,'d')), 16)

y1specific = 

The values near 1e-23 might arithmetically be 0. On your system it might be worth vpa to more digits than the 16 I use here, and might be worth using sym(value_of_constants) instead of sym(value_of_constants,'d')

Although there appears to be an exact solution in algebraic numbers and exp() of algebraic numbers, it would be an extremely long solution.

Answer 3

[Sam Chak]

Hi @Chris G.

If you're using the ode45 solver to find the state responses, here is how they look. Indeed, as mentioned by @Walter Roberson, the analytical solution for is very lengthy, and you can find it here on WolframAlpha.

tspan  = [0, 16];
y0     = [0; 0; 0];
[t, y] = ode45(@odefcn, tspan, y0);

% Solution Plot
plot(t, y, 'linewidth', 1.5), grid on, xlabel('t'), ylabel('\bf{y}')
legend('y_{1}', 'y_{2}', 'y_{3}', 'fontsize', 12)

function ydot = odefcn(t, y)
    ydot    = zeros(3, 1);
    f       =  10*exp(-t);
    fprime  = -10*exp(-t);
    ydot(1) = y(2);
    ydot(2) = y(3);
    ydot(3) = - 7*y(3) - 8*y(2) - 9*y(1) + fprime + 2*f + 3*fprime + 4*f;
end

[Sam Chak]

Relatively short analytical solutions for a 3rd-order LTI system exist when the eigenvalues of the system do not contain complex values.

% characteristic polynomial: y''' + 6·y'' + 11·y' + 6·y = 0
p     = [1 6 11 6]; 
r     = roots(p)    % eigenvalues

r = 3×1

-3.0000 -2.0000 -1.0000

syms t y1(t) y2(t) y3(t)

f     = 10*exp(-t);
eqn1  = diff(y1) == y2;
eqn2  = diff(y2) == y3;
eqn3  = diff(y3) == - 6*y3 - 11*y2 - 6*y1 + diff(f,3) + 2*diff(f,2) + 3*diff(f, 1) + 4*f;
eqns  = [eqn1, eqn2, eqn3];
cond1 = y1(0) == 0;
cond2 = y2(0) == 0;
cond3 = y3(0) == 0;
conds = [cond1, cond2, cond3];
y     = dsolve(eqns, conds, 'MaxDegree', 3)

y = struct with fields:
    y2: (exp(-2*t)*(80*exp(t) - 80))/2 - 10*t*exp(-t) - (exp(-3*t)*(45*exp(2*t) - 45))/3
    y1: 10*t*exp(-t) - (exp(-2*t)*(80*exp(t) - 80))/4 + (exp(-3*t)*(45*exp(2*t) - 45))/9
    y3: 10*t*exp(-t) - exp(-2*t)*(80*exp(t) - 80) + exp(-3*t)*(45*exp(2*t) - 45)