Out of memory - dsolve 3rd order system

Question

Chris G. 2023-10-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2040236-out-of-memory-dsolve-3rd-order-system

评论： Sam Chak 2023-10-30

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-10-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2040236-out-of-memory-dsolve-3rd-order-system#answer_1343036

编辑：Torsten 2023-10-29

在 MATLAB Online 中打开

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)

2 个评论
显示无隐藏无

Chris G. 2023-10-30

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 2023-10-30

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

Or a numerical solution.

请先登录，再进行评论。

Answer 2

Walter Roberson 2023-10-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2040236-out-of-memory-dsolve-3rd-order-system#answer_1343046

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Walter Roberson 2023-10-30

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 2023-10-30

在 MATLAB Online 中打开

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 2023-10-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2040236-out-of-memory-dsolve-3rd-order-system#answer_1343156

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sam Chak 2023-10-30

在 MATLAB Online 中打开

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)