Solution of ODE with Symbolic Math Toolbox.

Hello,

I want to solve the ordinary differential equation

with , , . The function and should be a unity step. In other words: it is the step response of a first order lag element. I am using the Symbolic Math Toolbox.

First I define the variables including the relevant assumptions.

syms("t", "positive")
syms("T_1", "k_p", "positive")
syms("u(t)", "y(t)")

The I establish the differential equation

deq = T_1 * diff(y, t) + y == k_p * u

deq(t) = 

and replace with the unist step function.

deq = subs(deq, u, heaviside(t))

deq(t) = 

The solution

y_sol_t = dsolve(deq, y(0) == 0)

y_sol_t = 

y_sol_t = simplify(y_sol_t)

y_sol_t = 

is not very clear for students. It is hard to see that for .

If I use the Laplace Transform instead I get: the following result.

Laplace Transform of the above differential equation.

syms("s")
aeq = laplace(deq, t, s)

aeq = 

syms("y_LT")

Replacing laplace(y(t), t, s) and the initial condition

aeq = subs(aeq, [laplace(y(t), t, s) y(0)], [y_LT, 0])

aeq = 

Solve the algbraic equation vor y_LT

y_LT = solve(aeq, y_LT)

y_LT = 

Inverse Laplace Transform

y_sol_LT = ilaplace(y_LT, s, t)

y_sol_LT = 

which is valid only for , but can be extended to the general solt ion by multiplying it with ethe unit step.

y_sol_LT = y_sol_LT * heaviside(t)

y_sol_LT = 

How can I obtain this solution from the result calculated using dsolve?

Michael

Answer 1

[Star Strider]

One approach simply tells simplify to keep going until it reaches the set limit (500 here), or can't simplify further.

Try something like this --

syms("t", "positive")
syms("T_1", "k_p", "positive")
syms("u(t)", "y(t)")

deq = T_1 * diff(y, t) + y == k_p * u

deq(t) = 

deq = subs(deq, u, heaviside(t))

deq(t) = 

y_sol_t = dsolve(deq, y(0) == 0)

y_sol_t = 

y_sol_t = simplify(y_sol_t, 500)

y_sol_t = 

syms("s")
aeq = laplace(deq, t, s)

aeq = 

syms("y_LT")

aeq = subs(aeq, [laplace(y(t), t, s) y(0)], [y_LT, 0])

aeq = 

aeq = aeq * laplace(heaviside(t))

aeq = 

y_LT = solve(aeq, y_LT)

y_LT = 

y_sol_LT = ilaplace(y_LT, s, t)

y_sol_LT = 

y_sol_LT = y_sol_LT * heaviside(t)

y_sol_LT = 

Proof = isAlways(y_sol_t == y_sol_LT)

Warning: Unable to prove '-(k_p*(sign(t) + 1)*(exp(-t/T_1) - 1))/2 == heaviside(t)*(k_p - k_p*exp(-t/T_1))'.

Proof = logical
   0

Plotting --

(Since isAlways cannot prove that they are the same, plotting is appropriate.)-

T_1 = rand

T_1 = 0.4337

k_p = rand

k_p = 0.6785

figure
fplot(subs(y_sol_t), [-10 10])
grid

figure
fplot(subs(y_sol_LT), [-10 10])
grid

figure
fplot(subs(y_sol_t), [-10 10], LineWidth=3)
hold on
fplot(subs(y_sol_LT), [-10 10], LineWidth=1)
hold off
grid

.

Answer 2

[Paul]

Use rewrite to get the solution in terms of heaviside.

syms("t", "positive")
syms("T_1", "k_p", "positive")
syms("u(t)", "y(t)")
deq = T_1 * diff(y, t) + y == k_p * u;
deq = subs(deq, u, heaviside(t));
y_sol_t = dsolve(deq, y(0) == 0)

y_sol_t = 

y_sol_t = simplify(y_sol_t)

y_sol_t = 

y_sol_t = rewrite(y_sol_t,'heaviside')  % makes clear the solution is 0 for t < 0

y_sol_t = 

To get the exact same result as from Laplace ... an easy, if not elegant, approach

simplify(y_sol_t/heaviside(t))*heaviside(t)

ans =