%% Question 9 Solution
%Forward Transforms
syms t
laplace(exp(-t))
laplace(sinh(t))
s = linspace(1,10,100);
f_hat_exp = numeric_laplace(@(x)exp(-x),s);
f_hat_sinh = numeric_laplace(@(x)sinh(x),s);
figure(1)
plot(s,1./(s+1))
hold on
plot(s,f_hat_exp,’.’)
hold off
figure(2)
plot(s,1./(s.^2-1))
hold on
plot(s,f_hat_sinh,’.’)
hold off
%% Inverse Transforms
syms s
ilaplace((s-1)./((s-1).^2+1)-1./s)
ilaplace(2./s-1./(s+1))
t = linspace(0,5,100);
f_1 = numeric_inverse_laplace(@(x)(x-1)./((x-1).^2+1)-1./x,t,15);
f_2 = numeric_inverse_laplace(@(x)2./x-1./(x+1),t,10);
figure(3)
plot(t,exp(t).*cos(t)-1)
hold on
plot(t,f_1,’.’)
hold off
figure(4)
plot(t,2-exp(-t))
hold on
plot(t,f_2,’.’)
hold off
%%
%
% Code
function f_hat = numeric_laplace(f,s)
f_hat = zeros(size(s));
for i = 1:length(s)
dt = 0.01/s(i);
f_hatt = (dt/2)*(f(0));
f_hat_old = 100;
tol = 10e-10;
k = 1;
while abs(f_hatt-f_hat_old)>tol
f_hat_old = f_hatt;
f_hatt = f_hatt + dt*exp(-s(i)*k*dt)*f(k*dt);
k = k+1;
end
f_hatt = f_hatt + dt*exp(-s(i)*k*dt)*f(k*dt);
f_hat(i) = f_hatt;
end
end
%%
%
function f = numeric_inverse_laplace(f_hat,t,A,tol)
if nargin==3
tol = 10^-6;
end
f = zeros(1,length(t));
for m = 1:length(t)
ft = (exp(A/2)/(2*t(m)))*real(f_hat(A/(2*t(m))));
k = 1;
f_old = 1000;
while abs(ft-f_old)>tol
f_old = ft;
ft = ft+(exp(A/2)/t(m))*(-1)^(k)*real(f_hat((A+2*k*pi*1i)/(2*t(m))));
k = k+1;
end
f(m) = ft;
end
end