solving ODE by analytical & numerical solution

I'm a newbie to matlab and im trying to solve an ODE analytically and numerically

The problem given is x''+0.5x'+x=0

clc; clear all;

syms x(t) t s

m=1; b=0.5; k=1;

t_final=2;

Dx = diff(x,t);

D2x = diff(x,t,2);

F1 = laplace(m*D2x+b*Dx+k*x,t,s)

x(t) = ilaplace((s+0.5)/(s^2+0.5*s+1)) %analytical solution

sys = tf([1 0.5],[1 0.5 1]) %numerical solution

First I tried to laplace transform the equation and since I dont know how to automatically inverse the eq. by code, I had to organize the eq. manually and just used ilaplace to inverse the eq. and get the soltion.

For the numerical solution, I actually had no idea how to solve the problem that way so I just copied the code that was on an example problem but failed to get the answer.

This was the example problem that I refered. It was a 'car speed cruise control model'. In this code, they just used sys=tf(1, [m b]) to get a numerical solution.

My script just plots the eq on the laplace plane and doesnt match with the analytical solution

So my firtst question is are there any way to get a numerical solution just done by code? How can I get the Laplace transformed eq. to get organized and inversed back to get the answer by the code itself.

My second question is why does that single line of code [ sys = tf() ] works on the example script but doesnt work on mine. Perhaps can I be having the wrong toolbox? I have the "control system toolbox" added on.

[Dyuman Joshi]

Is it necessary to use laplace transform to solve the ODE?

[승언 김]

The class that Im taking requires us to learn several ways of handling ODE problems.

One of them is using Lapace transform to solve ODEs analytically.

So for me, yes!

Answer 1

[Sam Chak]

Hi @승언 김

There are a few approaches to solving linear differential equations without the Control System Toolbox. Here is one of them. The ode object numerical method was introduced in release R2023b.

%% Analytical Method
syms x(t)
Dx      = diff(x,1);
eqn     = diff(x,2) + 0.5*Dx + x == 0;
cond    = [x(0)==1, Dx(0)==0];
xSol(t) = dsolve(eqn, cond)

xSol(t) = 

fplot(xSol, [0 20]), grid on, xlabel('t'), ylabel('x(t)')
title('Analytical Method using dsolve')

%% Numerical Method
F              = ode;
F.ODEFcn       = @(t, x) [x(2);         % x1' = x2
                    - 0.5*x(2) - x(1)]; % x2' = - 0.5*x1' - x1
F.InitialValue = [1 0];
sol = solve(F, 0, 20)

sol = 
  ODEResults with properties:

        Time: [0 5.0238e-05 1.0048e-04 1.5071e-04 2.0095e-04 4.5214e-04 7.0333e-04 9.5452e-04 0.0012 0.0025 0.0037 0.0050 0.0062 0.0125 0.0188 0.0251 0.0313 0.0627 0.0941 0.1255 … ] (1×125 double)
    Solution: [2×125 double]

plot(sol.Time, sol.Solution(1,:), "-o"), grid on, xlabel('t'), ylabel('x(t)')
title('Numerical Method using ODE object')

[Sam Chak]

Q1: My script just plots the eq on the laplace plane and doesn't match with the analytical solution. So my first question is, are there any way to get a numerical solution just done by code?

The reason for the discrepancy between the numerical solution and the analytical solution is due to the use of the 'step()' command for plotting the step response. The 'step()' command assumes that the system starts at rest with x(0) = 0, and not at x(0) = 1. If the initial condition is known to be different, you should use the 'initial()' command or the 'lsim()' command. In the case of a time-variant external input signal, the 'lsim()' command is the appropriate choice.

Control System Toolbox: For linear ordinary differential equations (aka linear systems)

https://www.mathworks.com/help/control/ref/dynamicsystem.initial.html

https://www.mathworks.com/help/control/ref/dynamicsystem.lsim.html

For general ordinary differential equations, where the highest-order derivatives are linear:

https://www.mathworks.com/help/matlab/ref/ode45.html

[Sam Chak]

Q2: My second question is why does that single line of code "sys = tf()" works on the example script but doesnt work on mine. Perhaps can I be having the wrong toolbox?

The reason your script didn't work is because the transfer function of the system was specified incorrectly. The simulation problem aims to generate the time response to an initial condition without any external input signal. The step response corresponds to the system's response to the external Heaviside step signal, . The Laplace transform of the Heaviside function is . Therefore, to use the 'step()' command to simulate an input-free system, you need to modify the numerator of the transfer function so that the 's' variable will cancel out in the 'step()' command later.

syms x(t) t s X(s)

% parameters
m    = 1; 
b    = 0.5; 
k    = 1;

%% Analytical solution using Inverse Laplace method
Dx   = diff(x,t);
D2x  = diff(x,t,2);
Eqn  = m*D2x + b*Dx + k*x == 0

Eqn(t) = 

LEqn = laplace(Eqn);
LEqn = subs(LEqn, {laplace(x(t), t, s), x(0), Dx(0)}, {X(s), 1, 0});
LEqn = isolate(LEqn, X(s))

LEqn = 

x(t) = ilaplace(rhs(LEqn))

x(t) = 

fplot(x, [0 20]), grid on, xlabel('t'), ylabel('x(t)'), ylim([-0.5 1])
title('Analytical Method using ilaplace()')

%% Numerical solution using Step Response
Glap   = tf([2 1], [2 1 2])

Glap =
 
     2 s + 1
  -------------
  2 s^2 + s + 2
 
Continuous-time transfer function.

s      = tf('s');
Gaug   = s*Glap     % Augmented transfer function

Gaug =
 
    2 s^2 + s
  -------------
  2 s^2 + s + 2
 
Continuous-time transfer function.

tFinal = 20;

step(Gaug, tFinal), grid on, xlabel('t'), ylabel('x(t)')
title('Response to Initial Condition using step()')

[승언 김]

Thank you so much. You are my savior!!!

If you dont mind can I ask you another question?

I have another problem to be done

0.1*x' + x = 1(t), x(0) = 0

I understood how the Analytical method part works.

syms x(t) t s
Dx      = diff(x,1);
eqn     = 0.1*Dx + x == 1;
cond    = [x(0)==0];
xSol(t) = dsolve(eqn, cond)

xSol(t) = 

F              = ode;
F.ODEFcn       = @(t, x) [x(1);         
                    - 10*x(1) + 10];     %I meant it to be x1' = -10*x1 + 10
F.InitialValue = [0];
sol = solve(F, 0, 20)

Error using odearguments
@(T,X)[X(1);-10*X(1)-1] returns a vector of length 2, but the length of initial conditions vector is 1. The vector returned by @(T,X)[X(1);-10*X(1)-1] and the initial conditions vector must have the same number of elements.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in matlab.ode.internal.LegacySolver/solve (line 169)
[t,y] = obj.SolverFcn(obj.ODEFcn,[obj.InitialTime,t2], ...

Error in ode/solveBetween (line 699)
[tR,yR] = solver.solve(t0,tmax,pps);

Error in ode/solve (line 370)
sol = obj.solveBetween(t1,t2,nv.Refine);

But I couln't get the answer using Numerical method.

I would really appreciate if you could explain more specifically about using "F.ODEFcn" so that later on I could use and apply this code to other problems.

Answer 2

[Sam Chak]

You are welcome, @승언 김. If you find the solutions helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Thanks a bunch!

By the way, the reason your script for second ODE problem didn't work is because the first-order ODE function was incorrectly specified as a second-order system, , . I have corrected it in 'F.ODEFcn'. Take a look below.

syms x(t) t s
Dx      = diff(x,1);
eqn     = 0.1*Dx + x == 1;
cond    = [x(0)==0];
xSol(t) = dsolve(eqn, cond)

xSol(t) = 

tspan   = [0 1];
fplot(xSol, tspan), grid on, xlabel('t'), ylabel('x(t)')
title('Analytical Method using dsolve')

F              = ode;
F.ODEFcn       = @(t, x) - 10*x + 10;   % x' = - 10*x + 10
F.InitialValue = [0];
tFinal         = 1;
sol = solve(F, 0, tFinal)

sol = 
  ODEResults with properties:

        Time: [0 5.0238e-06 1.0048e-05 1.5071e-05 2.0095e-05 4.5214e-05 7.0333e-05 9.5452e-05 1.2057e-04 2.4616e-04 3.7176e-04 4.9735e-04 6.2295e-04 0.0013 0.0019 0.0025 0.0031 … ] (1×65 double)
    Solution: [0 5.0236e-05 1.0047e-04 1.5070e-04 2.0093e-04 4.5204e-04 7.0308e-04 9.5406e-04 0.0012 0.0025 0.0037 0.0050 0.0062 0.0124 0.0186 0.0248 0.0309 0.0608 0.0898 0.1180 … ] (1×65 double)

plot(sol.Time, sol.Solution, "-o"), grid on, xlabel('t'), ylabel('x(t)')
title('Numerical Method using ODE object')

[승언 김]

I really appreciate it. Thanks!!!