1

I have an MATLAB script that solves an inhomogeneous first-order linear IVP using the Laplace transform. (For this example, the script is set up to solve the IVP enter image description here, enter image description here.)

syms x(t) s X;

a0 = -3;
x0 = 4;
rhs = t^2;

lhs = diff(x,t) + a0*x;
ode = lhs - rhs

Lx = X;
LDx = s*X - x0;
LHS = LDx + a0*Lx;
RHS = laplace(rhs,t,s);
IVP = LHS - RHS;

IVP = collect(IVP,X);

X = solve(IVP, X);
X = partfrac(X);

sol = ilaplace(X, s, t)
check1 = diff(sol,t) - 3*sol
check2 = vpa(subs(sol, t, 0))

If I substitute "factor" for "collect", the script almost works on Octave with the symbolic package linking to SymPy, except for the "solve" command https://www.mathworks.com/help/symbolic/solve.html.

Is there any Octave (or SymPy, if that would function as a workaround) command that will function as a MATLAB symbolic toolbox "solve" command so I can solve the IVP with a Laplace transform with a script so I don't have to solve for X manually, then use "ilaplace"?

Thanks in advance for any assistance you can provide.

Tasos Papastylianou
  • 21,371
  • 2
  • 28
  • 57
  • What do you mean by "almost" works? Btw, `collect` [seems](https://uk.mathworks.com/help/symbolic/collect.html) to be used to collect coefficients; I spotted a `coeffs` function in the symbolic package which seems to be equivalent (even though, yes, `factor` seems to work in this case). – Tasos Papastylianou May 05 '17 at 08:49
  • based on the error I got when I ran this code (`Python exception: AttributeError: MutableDenseMatrix has no attribute is_Relational`), this might be relevant: http://stackoverflow.com/questions/42802588/solving-two-non-linear-equations-in-octave – Tasos Papastylianou May 05 '17 at 09:02
  • Thanks so much for replying! The coeffs function does indeed appear to be in what I'm interested! A simple matrix multiplication by [1;X] reproduces MATLAB's collect. I have some half-a**ed scripts that reproduce the above MATLAB scripts (approximately, at least), but I'll post those in an "answer" below. Thanks so much again for replying! – Jeffrey Rolland May 06 '17 at 19:39

2 Answers2

1

Here are some Octave scripts that (approximately, at least) reproduce the above MATLAB scripts. You have to manually enter the solution to the IVP = 0 in terms of X into Part 2 of each script, but they do function. If anyone has anyway of having Octave solve IVP = 0 in terms of X as in the MATLAB solve function, I would be glad to hear it.

This pair solves $\dot{x} - 3x = t^2$, $x(0) = 4$.

Part 1:

syms x(t) s X;

a0 = -3;
x0 = 4;
rhs = t^2;

lhs = diff(x,t) + a0*x;
ode = lhs - rhs

Lx = X;
LDx = s*X - x0;
LHS = LDx + a0*Lx;
RHS = laplace(rhs,t,s); % The t and s in laplace aren't necessary, as they are default
IVP = LHS - RHS;

coeff = coeffs(IVP,X);
IVP = coeff*[1;X]

Part 2:

syms x(t) s X;

X = -1*((-4*s^3-2)/s^3)/(s-3)

X = partfrac(X);

sol = ilaplace(X, s, t)
check1 = diff(sol,t) - 3*sol
check2 = vpa(subs(sol, t, 0))

This pair solves $\ddot{x} - 2\dot{x} - 3x = t^2$, $x(0) = 4$, $\dot{x}(0) = 5$.

Part 1:

syms x(t) s X;

a1 =-2;
a0 = -3;
x0 = 4;
xdot0 = 5;
rhs = t^2;

Dx = diff(x,t);
D2x = diff(x,t,2);
lhs = D2x + a1*Dx + a0*x;
ode = lhs - rhs

Lx = X ;
LDx = s*X - x0;
LD2x = s^2*X - x0*s - xdot0;
LHS = LD2x + a1*LDx + a0*Lx;
RHS = laplace(rhs,t,s); % The t and s in laplace aren't necessary, as they are default
IVP = LHS - RHS;

coeff = coeffs(IVP,X);
IVP = coeff*[1;X]

Part 2:

syms x(t) s X;

a1 = -2;
a0 = -3;

X = -1*((-4*s^4 + 3*s^3 - 2)/s^3)/(s^2 - 2*s - 3)

X = partfrac(X);

sol = ilaplace(X, s, t)

Dsol = diff(sol,t);
D2sol = diff(sol,t,2);
check1 = D2sol + a1*Dsol + a0*sol
check2 = vpa(subs(sol, t, 0))
check3 = vpa(subs(Dsol, t, 0))

Thanks so much for all the help and suggestions! It is truly appreciated!

0

OK, so, one of my students solved the problem (I will contact him later in the week to see if he wishes to be publicly acknowledged regarding his solution).

You just need to define the result of coeff*[1;X] as an equation set equal to 0, say IVPEQ = coeff*[1;X] == 0, then use the symbolic package command solve on this equation, X = solve(IVPEQ, X).

Here is a version of my previous 1st-order IVP solver with my student's modification

syms x(t) s X;

a0 = -3;
x0 = 4;
rhs = t^2;

lhs = diff(x,t) + a0*x;
ode = lhs - rhs

Lx = X;
LDx = s*X - x0;
LHS = LDx + a0*Lx;
RHS = laplace(rhs,t,s); % The t and s in laplace aren't necessary, as they are default
IVP = LHS - RHS;

coeff = coeffs(IVP,X);
IVPEQ = coeff*[1;X] == 0;

X = solve(IVPEQ,X);

X = partfrac(X);

sol = ilaplace(X, s, t)

Dsol = diff(sol,t);
check1 = Dsol + a0*sol
check2 = vpa(subs(sol, t, 0))

and here is the 2nd-order IVP solver with the student's modification

syms x(t) s X;

a1 =-2;
a0 = -3;
x0 = 4;
xdot0 = 5;
rhs = t^2;

Dx = diff(x,t);
D2x = diff(x,t,2);
lhs = D2x + a1*Dx + a0*x;
ode = lhs - rhs

Lx = X ;
LDx = s*X - x0;
LD2x = s^2*X - x0*s - xdot0;
LHS = LD2x + a1*LDx + a0*Lx;
RHS = laplace(rhs,t,s); % The t and s in laplace aren't necessary, as they are default
IVP = LHS - RHS;

coeff = coeffs(IVP,X);
IVPEQ = coeff*[1;X] == 0;

X = solve(IVPEQ,X);

X = partfrac(X);

sol = ilaplace(X, s, t)

Dsol = diff(sol,t);
D2sol = diff(sol,t,2);
check1 = D2sol + a1*Dsol + a0*sol
check2 = vpa(subs(sol, t, 0))
check3 = vpa(subs(Dsol, t, 0))

Thanks again, @Tasos_Papastylianou , for your immense help!