4 LESSON 3
4.1 ORDINARY DIFFERENTIAL EQUATION
An Ordinary Differential Equation (ODE) is an equation (or equations system) that contains
one or more derivatives with respect to an independent variable (𝑡). MATLAB contains several
solvers and methods for the analytical, numerical and discrete solution of the ODE [15], [16].
The numerical solvers can solve stiff or nonstiff problems, problems with a mass matrix,
differential algebraic equations (DAEs), or fully implicit problems. The numerical solver
selection for ODE is describe in [15].
We are looking for a solution of simple an ODE which describes an RL series electrical circuit
(Fig. 4-1). The solution of an ODE is the current waveform, depending on the initial conditions.
i(t)
u(t) R L
Fig. 4-1. RL series electric circuit: u – voltage source, R –resistance, L – inductance, i – current.
The voltage equation of an RL series electric circuit is defined as:
d𝑖(𝑡)
𝑢(𝑡) = 𝑅 ⋅ 𝑖(𝑡) + 𝐿 ⋅
d𝑡
Where: 𝑢(𝑡) – supply voltage, 𝑅 – resistance, 𝐿 – inductance, 𝑖(𝑡) – current.
4.1.1 Symbolic solution
Symbolic Math Toolbox from MATLAB provides functions for solving, plotting, and
manipulating symbolic math equations. The equivalent of this package for Octave is The
Octave-Forge Symbolic package [17]. Both of them can be used to find a symbolic solution of
an ODE.
To find a symbolic solution of an RL series equation, first, the symbolic variables must be
defined. To create symbolic variables use syms function [16]:
syms i(t) R L U; % symbolic variables
To define the differential equation use double equality symbol ‘==’ and diff as the derivative
symbol.
eqn = L*diff(i,t)+R*i == U;
Solve the ODE equation using the dsolve function.
26/45
tic % Start stopwatch timer
y=dsolve (L*diff(i) + R*i == U) % Define differential equation ‘diff’
toc % Stop stopwatch timer
pretty(y); % Prettyprint symbolic expressions
Rusults:
y =
(U - C2*exp(-(R*t)/L))/R
Elapsed time is 0.049207 seconds.
/ R t \
U - C2 exp| - --- |
\ L /
-------------------
R
The constant C2 appears because no initial condition was specified. You can find the value of
constant by solve the equation with the initial condition.
eqn = L*diff(i,t)+R*i == U
ic = 'i(0)=U/R'
y=dsolve(eqn, ic) % initial condition i(0)=U/R
y =
(U - U*exp(-(R*t)/L))/R
The result of symbolic calculations can also be presented in the form of a plot. For this purpose,
the symbolic variable should be assigned parameter values and the range of the independent
variable (time) should be defined.
Check done for ODE solution with initial condition and i(0) = 0.
clc;
clear;
close all;
%%% ODE solution
syms i(t) R L U;
eqn = L*diff(i,t)+R*i == U
y=dsolve(eqn, 'i(0)=0')
%%% Input data
U=10; % [V] supply voltage
R=1; % [Ohm] resistance
27/45
L=1e-3; % [H] inductance
T=L/R; % [s] time constant
t=0:T/100:5*T; % [s] time range
i= subs(y); % Symbolic substitution – each occurrence of the variable i
%
%%%
% Displaying the graph with the ODE solution
%
figure; % graphic window
plot(t,i,'r-'); % plot results
xlabel('time t[s]'); % x-axis label
ylabel('current i[A]'); % y-axis label
title('Series RL'); % plot title
grid on; % grid on
Fig. 4-2. Resutls of analytical solution of the ODE
28/45
4.1.2 Numerical solution
To find a numerical solution of the ODE appropriate solvers should be used. The solver Runge-
Kutta (4,5) name ode45 can be used to solve most ODEs. More efficient and faster solvers to
use is ode23 (less accurate) and ode113 (bigger problems, tighter accuracy requirements).
All numerical MATLAB ODE solvers can solve systems of equations of the form:
𝑦 ′ = 𝑓(𝑡, 𝑦)
or problems that involve a mass matrix:
𝑀(𝑡, 𝑦)𝑦 ′ = 𝑓(𝑡, 𝑦)
All ODE solvers use similar syntaxes.
[t,y] = solver (@fun, tspan, init,option)
Where:
• [t, y] vectors with calculation results,
• @fun or 'fun' - a function handle or function name containing the ODE definition
• tspan – time vector,
• init – initial conditions,
• option – solver parameters and other options.
clc; % clear command window
clear; % clear workspace memory
close all; % close all graphic windows
%%% Input data
U=10; % [V] supply voltage
R=1; % [Ohm] resistance
L=1e-3; % [H] inductance
T=L/R; % [s] time constant
tspan=0:T/100:5*T; % [s] time range
i0=[0]; % [A] initial condition
%%% ODE solution
tic
[t1,i1] = ode45('series_RL',tspan,i0);
toc
As first argument for the ode45 function the a function handle or function name containing the
ODE definition is used. In this example, the new file with function ‘series_RL’ should be
created.
function didt = series_RL(t,i)
%SERIES_RL
U=10; % [V] supply voltage
R=1; % [Ohm] resistance
29/45
L=1e-3; % [H] inductance
%%%
% ODE equation:
didt=zeros(1);
didt = 1/L*U - R/L*i(1);
end
Passing the name of ODE function in this way is inconvenient because you can not pass other
parameters, eg values of model parameters. A more flexible approach is provided by the
definition of an inline function.
% inline function for passing parameters
series_RL=@(t,i)series_RL_param(t,i,U,R,L);
tic
[t2,i2] = ode45(series_RL,tspan,i0);
toc
function didt = series_RL_param(t,i,U,R,L)
%%%SERIES_RL_param
%
didt=zeros(1); % Resetting result vector / matrix
didt = 1/L*U - R/L*i(1); % derivative of the function
end
% Visualization of the ODE solution
figure; % new graphic window
plot(t1,i1,'r-'); % plot first solution
hold on; % hold on
plot(t2,i2,'bo'); % plot second solution
xlabel('czas t[s]'); % X-axis label
ylabel('prąd i[A]'); % Y-axis label
title('RL series circuit'); % title
legend('solution i_1','solution i_2');
grid on;
30/45
Fig. 4-3. Resutls of numerical solution of the ODE
EX.1
Find a solution for an RL series electric circuit with the sinusoidal voltage
supply.
4.1.3 Discrete data
clc; % clear command window
clear; % clear workspace memory
close all; % close all graphic windows
% Generate discrete data
x=0:pi/10:2*pi; % x range
y=sin(x); % discrete function
dydx=diff(y)./diff(x); % Derivative of discrete data
Ny = length(y); % Vector length y => N
Ndy = length(dydx); % Vector length dydx => N-1
disp([Vector length Ny=', num2str(Ny)]);
disp([Vector length Ndy=', num2str(Ndy)]);
% Resutls
figure; % Figure window divided into two graphs
31/45
subplot(2,1,1); % subplot no 1
plot(x,y,'rx-'); % plot
xlim([0 2*pi]); % Set x-axis limits
xlabel('x'); % x-axis label
ylabel('y=f(x)'); % y-axis label
grid on; % grid switch on
subplot(2,1,2); % subplot 2
plot(x(1:end-1),dydx,'mo-');
xlim([0 2*pi]);
xlabel('x');
ylabel('dy/dx=f''(x)');
grid on;
4.2 INTEGRATION
4.2.1 Symbolic solution
clc; % clear command window
clear; % clear workspace memory
close all; % close all graphics windows
%% Symbolic integration
syms x % symbolic variable
y = int(sin(x)) % indefinite integral
pretty(y); % pretty solution
Y0=int(sin(x),0,pi) % integral define in range 0..pi
4.2.2 Numerical solution
Integral numerically integrates function fun from xmin to xmax using global adaptive quadrature and
default error tolerances
Trapz - Trapezoidal numerical integration
clc; % clear command window
clear; % clear workspace memory
close all; % close all graphics windows
x=0:pi/10:2*pi; % x range
fun = @(x)sin(x); % inline function
Y1 = quad(fun,0,pi) % Simpson method
32/45
Y2 = integral(fun,0,pi) % Integral
%%%
% Porównaj wyniki
Y3 = trapz(x, sin(x)) % Trapezoidal numerical integration
Y4 = trapz(x, abs(sin(x))/2) % Trapezoidal numerical integration
x=0:pi/100:pi; % x range
Y5 = trapz(x, sin(x)) % Trapezoidal numerical integration
4.3 RLC CIRCUIT
The second order ODE will be analyzed using the example of a series RLC electric circuit. Fig.
4-4 shows the series form the RLC electric circuit.
i(t) R
u(t) L
C
Fig. 4-4. RLC series electric circuit
The voltage equations using :
𝑑
𝑢(𝑡) = 𝑅 ⋅ 𝑖(𝑡) + 𝐿 𝑖(𝑡) + 𝑢𝐶 (𝑡)
𝑑𝑡
Where
1
𝑢𝑐(𝑡) = ∫ 𝑖(𝑡)𝑑𝑡
𝐶
using substitution:
𝑑
𝑖(𝑡) = 𝑞(𝑡)
𝑑𝑡
the homogeneous second order ODE for the voltage equation is given by:
𝑑 𝑑2 1
𝑢(𝑡) = 𝑅 ⋅ 𝑞(𝑡) + 𝐿 ⋅ 2 𝑞(𝑡) + ⋅ 𝑞(𝑡)
𝑑𝑡 𝑑𝑡 𝐶
4.3.1 Symbolic solution
The symbolic solution of second order ODE with the initial condition.
33/45
Clear workspace memory
clc; % clear command window
clear; % clear workspace memory
close all; % close all graphics windows
Define the symbolic variables of ODE.
tic % timer start
syms q(t) R L C U; % symbolic variable
General form of second order ODE solution (without initial condition)
y=dsolve('L*D2q + R*Dq + 1/C*q ==U')
pretty(y);
ODE solution with initial conditions
Dq = diff(q);
y=dsolve(L*diff(q,2) + R*diff(q,1) + 1/C*q ==U, q(0)==0,Dq(0)==0);
Results
disp('Charge equation q(t)');
pretty(y); %
disp('Current equation i(t)=dq(t)/dt');
pretty(diff(y)); %
toc % timer stop
Model parameters
U=10; % [V] supply voltage
R=1e3; % [Ohm] resistance
L=1; % [H] inductance
C=0.5e-6; % [F] capacitance
Time constants
T_RL=L/R % [s] RL time constant
T_RC=R*C % [s] RC time constant
t=0:T_RL/100:10*T_RL; % [s] t range
Paramterers substitution
34/45
q = subs(y); % charge waveform q(t)
i = subs(diff(y)); % current waveform i(t
Graphical solution
figure; % new graphic figure
subplot(3,1,1);
plot(t,q,'r-' ); % charge waveform q(t)
title('charge q(t)' ); % plot title
xlabel('t[s]' ); % x-axis label
ylabel('q[C]' ); % y-axis label
grid on ; % grid
subplot(3,1,2);
plot(t,q/C,'g-' ); % capacitor voltage uc(t)
title('voltage uc(t)' ); % plot title
xlabel('t[s]' ); % x-axis label
ylabel('uc[V]' ); % y-axis label
grid on ; % grid on
subplot(3,1,3);
plot(t,i,'b-' ); % current waveform i(t)
title('current i(t)' ); %
xlabel('t[s]' ); %
ylabel('i[A]' ); %
grid on ;
Fig. 4-5. Results of second order ODE
35/45
4.3.2 Numerical solution
Numerical solution of second order ODE.
function dxdt = series_RLC_param(t,x,U,R,L,C)
dxdt=zeros(2,1);
dxdt(1) = 1/C*x(2);
dxdt(2) = 1/L*U - R/L * x(2) -1/L*x(1);
end
Clear workspace memory
clc; % clear command window
clear; % clear workspace memory
close all; % close all graphics windows
Model parameters
U=10; % [V] supply voltage
R=1e3; % [Ohm] resistance
L=1; % [H] inductance
C=0.5e-6; % [F] capacitance
Time constants
T_RL=L/R % [s] RL time constant
T_RC=R*C % [s] RC time constant
t=0:T_RL/100:10*T_RL; % [s] t range
Simulaiton parameters
• 'RelTol' — Relative error tolerance – This tolerance measures the error relative to the
magnitude of each solution component. Roughly speaking, it controls the number of
correct digits in all solution components, except those smaller than the absolute tolerance
AbsTol.
• 'AbsTol' — Absolute error tolerance - This tolerance is a threshold below which the value
of the solution becomes unimportant.
• 'NormControl' — Control error relative to norm
• 'OutputFcn' — Output function – the ODE solver calls the output function after each
successful time step.
x0 = [0,0]; % initial condition/values
tspan = [0 20e-3];
option = odeset ('RelTol', 1e-6, 'AbsTol', 1e-6, 'NormControl', 'on', ...
'OutputFcn', @odeplot);
series_RLC = @(t,x)series_RLC_param(t,x,U,R,L,C);
[t1,x1] = ode45(series_RLC,tspan,x0,option);
[t2,x2] = ode23(series_RLC,tspan,x0,option);
36/45
Fig. 4-6. Numerical solution control - Plot components of the solution vs. time
Plot results
figure;
subplot(2,1,1);
plot(t1,x1(:,1));
xlabel('time t[s]' );
ylabel('voltage vc[V]' );
grid on ;
subplot(2,1,2);
plot(t2,x2(:,2));
xlabel('time t[s]' );
ylabel('current i[A]' );
grid on ;
37/45
Fig. 4-7. Results of numerical solution of second order ODE
EX. 2
Analyze the solution for different circuit parameters
Analysis of a RLC model parameters
alpha=1/2*T_RL % constant damping
wn = 1/sqrt(T_RL*T_RC) % the natural frequency of undamped vibrations
D = alpha/wn % damping factor
w = wn*sqrt(1-D^2) % the natural frequency of damped vibrations
alpha = 5.0000e-04
wn = 1.4142e+03
D = 3.5355e-07
w = 1.4142e+03
Possible solutions:
• 𝐷 > 1 𝑇𝑅𝐶 > 4 ⋅ 𝑇𝑅𝐿 – aperiodic waveform
• 𝐷 = 1 𝑇𝑅𝐶 = 4 ⋅ 𝑇𝑅𝐿 – an aperiodic critical waveform
• 𝐷 < 1 𝑇𝑅𝐶 < 4 ⋅ 𝑇𝑅𝐿 – oscillation waveform
38/45
