Solve Complex Maths Like ODE with Laplace, Maclaurin's Series Expansion, Fourier Series Expansion and the Heat Equation (PDE) with Matlab. [Cardis!]
1. Solve d'' + 6 d' + 9y using Laplace
clc; clear; syms t s Y y(t) Dy(t)
assume([t Y] > 0)
Dy = diff(y, t)
D2y = diff(Dy, t)
LS = sin(3*t)
% Diff. Equation formulation:
EQN=D2y+6 * Dy+9 * y-LS
% Laplace transform of the Diff. Equation
LEQN=laplace(EQN,t,s)
% Substitute ICs and initiate the arbitrary unknown
"Y"
LT_Y=subs(LEQN,laplace(y,t,s),Y)
LT_Y=subs(LT_Y, y(0), 0) % y(0) = 1
LT_Y=subs(LT_Y, subs(diff(y(t), t), t, 0), 0) %
dy(0)= 0
% Solve for the arbitrary unknown: Y
ys=solve(LT_Y,Y)
% Inverse of the Laplace Transform:
y=ilaplace(ys,s,t)
Description:
%First we cleared the workspace, the variables and declared
variables as symbolic. %Declaring variables as symbolic is necessary in using
the symbolic toolbox where %we have most of the tools used in this exercise.
clc; clear; syms t s Y y(t) Dy(t)
%Next we assume initial value of t = 0
assume([t] == 0)
%Then, we declare y to be a differential term in respect to t
Dy = diff(y, t)
D2y = diff(Dy, t)
%Those two were rewritten as:
Output:
Dy(t) =
D2y(t) =
%… in the output.
%We declared the left hand side to be sin(t)
LS = sin(3*t)
%Which was rewritten as:
Output:
LS =
%… in the output.
% Differential Equation was then formulated with these
variables and its Laplace % transform was found
EQN=D2y+6 * Dy+9 * y-LS
LEQN=laplace(EQN,t,s)
Output:
LEQN =
% However, because Matlab only solves for symbolic variables, we
need to %substitute complex terms with symbolic ones using the sub function. We
also %substituted for our initial values of x = 0 and xI = 0 into
the resultant %equation.
% Substitute ICs and initiate the arbitrary unknown
"Y"
LT_Y=subs(LEQN,laplace(y,t,s),Y)
LT_Y=subs(LT_Y, y(0), 0) % y(0) = 1
LT_Y=subs(LT_Y, subs(diff(y(t), t), t, 0), 0) %
dy(0)= 0
Output:
LT_Y =
LT_Y =
LT_Y =
% Solve for the arbitrary unknown: Y
ys=solve(LT_Y,Y)
% Inverse of the Laplace Transform:
y=ilaplace(ys,s,t)
Output:
ys =
y =
% Final
function.
1.a.ii. There are many methods that can be used both
cursively and digitally. For instance, cursively, we can use the following:
i.
Method of undetermined coefficient
ii.
Integrating factors
iii.
Exact differential equation depending on the
circumstance.
Digitally, in Matlab, for instance, we can:
i.
Convert the ODE to a lower order ODE, and then
solve with the Laplace or ODE command
ii.
We may decide to solve directly with the dsolve
function, inputting the initial values as parameters in the command.
More so, there are many reasons to choose the laplace
technique when solving ODEs. Some of them are:
i.
It is relatively easier and faster. using
Laplace transforms reduces a differential equation down to an algebra problem.
This algebra, will often be easier than a straight forward approach. First, we
do not need to find a general solution, differentiate this, plug in the initial
conditions and then solve for the constants to get the solution. With Laplace
transforms, the initial conditions are applied during the first step and at the
end we get the actual solution instead of a general solution.
Furthermore, many equivalent values are already available in ready-made
tables which are now plug-and-play. For instance,
Laplace transform of dy/dx = sY(s) – Y(0) and so on…
A table example can be found in: https://tutorial.math.lamar.edu/classes/de/Laplace_Table.aspx
ii.
Finally, there are some differential equations
that simply can’t be done using the techniques from the last chapter and so, in
those cases, Laplace transforms will be our only solution.
2. Maclaurin’s series expansion of 
with percentage error for the first four non-zero terms.
syms x f
f = exp(x) * sin(x)
deriv0 = vpa(subs(f,x,0))
dfdx = diff(f)
deriv1 = vpa(subs(dfdx,x,0))
df2dx = diff(dfdx)
deriv2 = vpa(subs(df2dx,x,0))
df3dx = diff(df2dx)
deriv3 = vpa(subs(df3dx,x,0))
df4dx = diff(df3dx)
deriv4 = vpa(subs(df4dx,x,0))
df5dx = diff(df4dx)
deriv5 = vpa(subs(df5dx,x,0))
df6dx = diff(df5dx)
deriv6 =
vpa(subs(df6dx,x,0))
deriv = [deriv0 deriv1 deriv2 deriv3 deriv4 deriv5 deriv6]
t = sym(zeros());
for i = 0 : length(deriv)-1
% Implement the Maclaurin
expansion
t(i+1, :) = deriv(1) *
x.^(i) / factorial(i);
% Find the derivative
for the next term
deriv = circshift(deriv,
-1);
end
% Add-up the calculated terms
smp =
sum(t)
approx =
vpa(subs(smp,x,0.5))
exact =
vpa(subs(f,x,0.5))
p_error =
(abs(approx - exact)/abs(exact)) * 100
Description:
Both dependent and independent variables are declared as symbolic:
syms x f
Initial equation was then declared, after which all derivatives
were found for the four non-zero terms. Manual monitoring was done through the
outputs to ensure non-zero values:
f = exp(x) * sin(x)
deriv0 = vpa(subs(f,x,0))
dfdx = diff(f)
deriv1 = vpa(subs(dfdx,x,0))
df2dx = diff(dfdx)
deriv2 = vpa(subs(df2dx,x,0))
df3dx = diff(df2dx)
deriv3 = vpa(subs(df3dx,x,0))
df4dx = diff(df3dx)
deriv4 = vpa(subs(df4dx,x,0))
df5dx = diff(df4dx)
deriv5 = vpa(subs(df5dx,x,0))
df6dx = diff(df5dx)
deriv6 =
vpa(subs(df6dx,x,0))
Output:
f =
deriv0 =
dfdx =
deriv1 =
df2dx =
deriv2 =
df3dx =
deriv3 =
df4dx =
deriv4 =
df5dx =
deriv5 =
df6dx =
deriv6 =
All derivatives are stored in a vector for maclaurin’s expansion.
deriv = [deriv0 deriv1 deriv2 deriv3 deriv4 deriv5 deriv6]
t = sym(zeros());
for i = 0 : length(deriv)-1
% Implement the
Maclaurin expansion
t(i+1, :) = deriv(1) *
x.^(i) / factorial(i);
% Find the derivative
for the next term
deriv = circshift(deriv,
-1);
end
% Add-up the calculated terms
smp =
sum(t)
Output:
deriv =
smp =
Approximated
and exact values was computed for x = 0.5.
Finally,
percentage error was computed:
approx =
vpa(subs(smp,x,0.5))
exact =
vpa(subs(f,x,0.5))
p_error =
(abs(approx - exact)/abs(exact)) * 100
Output:
approx =
exact =
p_error =
1.b.ii.
Fourier series expansion
syms t w x n
%Fourier
series expansion of a periodic function is given by
%F(x) =
Ao + Eps(An * cos(n * x)) + Eps(Bn * sin(n * x))
%Where Ao
is the coefficient of the first term
%An is
the coefficient of the second term
%Bn is
the the coefficient of the third term
%Declaring
f(x) and their limits
f1 = -1
f1_lims = [-pi 0]
f2 = 1;
f2_lims = [0 pi]
% Compute first term using:
% Ao = 1/2pi * {int0-pi(f1)dx + int0pi
(f2)dx}
Ao = (1/(2*pi)) * (int(f1, x, f1_lims(1), f1_lims(2)) + int(f2, x,
f2_lims(1), f2_lims(2)))
% Compute second term:
% An = 1/pi * {int0-pi(f1)cosnxdx +
int0pi(f2)cosnxdx}
An = (1/pi) * (int(f1 * cos(n * x), x, f1_lims(1), f1_lims(2)) +
int(f2 * cos(n * x), x, f2_lims(1), f2_lims(2)))
% Compute third term:
% An = 1/pi * {int0-pi(f1)sinxndx +
int0pi(f2)sinxndx}
Bn = (1/pi) * (int(f1 * sin(n * x), x, f1_lims(1), f1_lims(2)) +
int(f2 * sin(n * x), x, f2_lims(1), f2_lims(2)))
% Add them all three terms using:
% F(x) =
Ao + Eps(An * cos(n * x)) + Eps(Bn * sin(n * x))
Fx = Ao + symsum(An .* cos(n .* x)) + symsum(Bn .* sin(n .* x))
Output:
f1
= -1
f1_lims
= 1×2 -3.1416 0
f2_lims
= 1×2 0
3.1416
Ao = 0
An = 0
Bn =
Fx =
2. The Heat
Equation
To
solve this equation in MATLAB, we coded the equation, initial conditions, and
boundary conditions, then select a suitable solution mesh before calling the
solver pdepe.
Code
Initial Condition
The
initial condition function for the heat equation assigns a constant value for .
This function must accept an input for x, even if it is unused.
function u0
= heatic(x)
u0 = 0.5;
end
Code
Boundary Conditions
function [pl,ql,pr,qr]
= heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end
Select
Solution Mesh
Use a
spatial mesh of 20 points and a time mesh of 30 points. Since the solution
rapidly reaches a steady state, the time points near are more closely spaced
together to capture this behavior in the output.
L
= 1;
x
= linspace(0,L,20);
t
= [linspace(0,0.05,20), linspace(0.5,5,10)];
Solve
Equation
Finally,
solve the equation using the symmetry m, the PDE equation, the
initial condition, the boundary conditions, and the meshes for x and
t.
m
= 0;
sol
= pdepe(m,@heatpde,@heatic,@heatbc,x,t)
Local
Functions
function
[c,f,s] = heatpde(x,t,u,dudx)
c
= 1;
f
= dudx;
s
= 0;
end
function
u0 = heatic(x)
u0
= 0.5;
end
function
[pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl
= ul;
ql
= 0;
pr
= ur - 1;
qr
= 0;
end
Comments
Post a Comment