Master Matlab

1.  Basics

Example: If A= [3 0 1, 0 1 0, 2 0 0] find (a) Characteristic polynomial of A, (b) Roots of the characteristic polynomial A, (c) eigenvalues of A, (d) eigenvectors of A, (e) find eigenvalues of A -1 , (f) find eigenvalues of A T , (g) find eigenvalues of B=A 2 +3A+2I
A=[3 0 -1; 0 1 0; 2 0 0]
Eigenvalues_of_A=eig(A)
[X,D]= eig(A)
p=poly(A) % Coefficients of the characteristic polynomial
r = roots(p) % To find the roots of the characteristic equation
sum_of_eigenvalues=sum(r)
Trace_of_A=trace(A) %Note that sum of eigenvalues is Trace of A
product_of_eigenvalues_A=prod(r)
Determinant_of_A =det(A) %Note that product of the eigenvalues of A is determinant of A
Inverse_A=inv(A)
Eigenvalues_of_invA= eig(Inverse_A)
Eigenvalues_Transpose_of_A= eig(A') % Note that eigenvalues of A and Transpose of A are same
Eigenvalues_of_B= eig(A^2+3*A+2*eye(3)) % Understand the property it is satisfying


Aim : Matrix Manipulation in Matlab

Q1
a=[1 2 3;4 5 6;7 8 9]
m=size(a)
a(3,2)

Q2
b=[1 1 3;1 5 1;3 1 1]
det(b)
rank(b)
trace(b)
diag(b)

Q3
b=[1 1 3;1 5 1;3 1 1]
b'
b^2
inv(b)
b^-1

Q4
c=a+b
d=a-b
e=a*b
f=a.*b
g=a./b

Q5
x=[a b]
y=[a;b]

A5

x =

     1     2     3     1     1     3
     4     5     6     1     5     1
     7     8     9     3     1     1


y =

     1     2     3
     4     5     6
     7     8     9
     1     1     3
     1     5     1
     3     1     1


Q6
poly(b)
roots(poly(b))
eig(b)
[v d]=eig(b)

A6

ans =
    1.0000   -7.0000   -0.0000   36.0000

ans =
    6.0000
    3.0000
   -2.0000

ans =
   -2.0000
    3.0000
    6.0000

v =
   -0.7071    0.5774    0.4082
         0   -0.5774    0.8165
    0.7071    0.5774    0.4082


d =
   -2.0000         0         0
         0    3.0000         0
         0         0    6.0000

Q7
x=[1 -1 0 1 2 3 4 5]
y=[0 1 2 3 4 5 6 7]
x.*y
max(x)
min(y)

Q8
a=1:5
a=0:2:10
a=linspace(1,5)
a=linspace(1,5,10)

Q9
ones(3)
ones(2,3)
zeros(3,3)
eye(2)

Q10
a=[4 -2 6;2 8 2;6 10 3]
b=[8;4;0]
inv(a)*b

Q11
a=[3 0 -1;0 1 0;2 0 0]
poly(a)
eig(a)
[v d]=eig(a)
eig(a^-1)
eig(a')
b=a^2+3*a+2*eye(3)
eig(b)


2.  Graph

Axis
>> x=[0:0.01:10];
>> y=exp(-x).*sin(2*x + 3);
>> plot(x,y),axis([0 10 -1 1])

Title and Label -
>> x=linspace(0,2*pi,100);
>> y=sin(x);
>> plot(x,y);
>> title('this is sinus function');
>> xlabel('x-values');
>> ylabel('y-values');

Multi Plot
>> x=[0:0.01:10];
>> y=sin(x);
>> g=cos(x);
>> plot(x,y,x,g,'.-'), legend('sin(x)','cos(x)')

Grid
>> t=0:pi/200:2*pi;
>> y1=sin(t);
>> y2=sin(t+pi/2);
>> plot(t,y1,t,y2);
>> grid on

Subplot 1
>> t=0:pi/200:2*pi;
>> y1=sin(t);
>> y2=sin(t+pi/2);
>> subplot(2,2,1)
>> plot(t,y1)
>> subplot(2,2,2)
>> plot(t,y2)

Subplot 2
>> x=[0:0.1:10];
>> f=figure;
>> f1=subplot(1,2,1);
>> plot(x,cos(x),'r');
>> grid on;
>> title('cosine');
>> f2=subplot(1,2,2);
>> plot(x,sin(x),'d');
>> grid on;
>> title('sine');

Polynomial Function
>> s=linspace(-5,5,100);
>> coeff = [1 3 3 1];
>> A = polyval(coeff,s);
>> plot(s,A)
>> title('Polynomial Function');
>> grid on;

3.  Caley Hamilton

clear all
clc
A=input('Enter the matrix A:');
n=length(A);
c=poly(A) % Coefficients of the characteristic polynomial
P=zeros(n,n); % Initialization purpose
for i=1:n+1
 P=P+c(i)*A^(n+1-i) % Characteristic equation
end
Determinant=det(A)
if det(A)~= 0
inv_A=zeros(n,n)
for j=1:n
 inv_A= inv_A - ((1/c(n+1))*(c(j)*A^(n-j)))
end
else
 disp( 'A is a singular matrix')
end

Input : A:[3 0 -1; 0 1 0; 2 0 0]


4.  Diagonalization
clear all;
clc;
A=input('Enter the matrix A:');
n=length(A);
[X,D]=eig(A);
fprintf('\n%s','The eigen values of A are')
disp(diag(real(D))')
disp('The eigen vectors for the corresponding eigenvalues')
disp(X)
option=input('If you want to diagonalize by similarity transformation then
press 1 or press any number:');
if(option==1)
 %Diagonalization of the real matrix A having distinct eigenvalues by
 %similarity transformation
 disp('Modal matrix associated with A is')
 P=X
 disp('D=inv(P)*A*P')
 D=inv(P)*A*P
 disp('Thus A is reduced to the diagonal matrix D through P by similarity
transformation.')
else
% Diagonalization of the real symmetric matrix A having distinct eigenvalues
by orthogonal transformation
for i=1:n,
 x=X(:,i);
 u(:,i)=x/norm(x);
end
disp('Orthogonal matrix associated with A is the matrix')
P1=u(:,1:n);
disp('D=transpose(P1)*A*P1')
D=(P1)'*A*(P1)
disp('Thus A is reduced to the diagonal matrix D through P1 by orthogonal
transformation.')
end


Output:
Enter the matrix A:[-1 2 -2;1 2 1;-1 -1 0]


5.  SOLUTION OF SYSTEM OF SECOND ORDER DIFFERENTIAL EQUATIONS OF THE FORM Y + AY = 0

clear all
clc
syms t a b
A = [1 2;2 -1];
l=eig(A);
[p,q]=eig(A);
s1 = dsolve('D2x+a*x=0','t')
s2 = dsolve('D2x+b*x=0','t')
disp(s1)
disp(s2)
disp('The Solution of the given system');
y = p*[s1;s2]
soln = subs(y,{a,b},{q(1),q(4)})

6.  SERIES SOLUTIONS

clear all
clc
syms x p n
n = input('Enter the number of terms up to which you want to evaluate the series (between 1 and 10) : ')
n=n-1;
disp('Note : Input the coefficients of DE as constant function of x');
p(1)= input('The coefficient of D2y : ');
p(2)= input('The coefficient of Dy : ');
p(3)= input('The coefficient of y : ');

for i = 0:n
    tmp = strcat('a', num2str(i));
    a(i+1) = sym(tmp);
end

y=sum(a.*x.^(0:n));
ps = collect(p(1)*diff(y,2) + p(2)*diff(y) + p(3)*y ,x);
psc=coeffs(ps,x);
i1=input('Enter the first IV: ');
i2=input('Enter the second IV: ');
ics1=strcat('a0=',num2str(i1));
ics2=strcat('a1=',num2str(i2));
sol=vpa(subs(y,solve(ics1,ics2,psc(1:(n-1)))),3);
disp('The solution of the given DE is:');
disp(sol);
ezplot(sol,[-4,4]);

7.  Bessels

Aim – Bessels Equation of first kind

X = 0:0.1:20;
J = zeros(5,nume1(X));
for i=0 : 4
    J(i+1,:) = besselj(i,X);
end
plot (X,J,'LineWidth',1.5)
grid on;
legend('J_0','J_1','J_2','J_3','J_4','Location','Best');
title('Bessel Functions of the First Kind for n=0,1,2,3,4')
xlabel('x')
ylabel('J_n(X)')



Aim – Bessels Equation of second kind

clc
clear all
X = 0:0.1:20;
Y = zeros(5,numel(X));
for i=0 : 4
    Y(i+1,:) = bessely(i,X);
end
plot (X,Y,'LineWidth',1.5)
axis([-0.1 20.2 -2 0.6]);
grid on;
legend('Y_0','Y_1','Y_2','Y_3','Y_4','Location','Best');
title('Bessel Functions of the First Kind for n=0,1,2,3,4')
xlabel('X')
ylabel('Y_n(X)')
8.  Fourier

clear all
clc
syms x k L U n
f = input('Enter the functions: ');
L = input('enter the lower limit: ');
U = input('Enter the upper limit: ');
l = (U-L)/2;
n = input('Enter the no. of terms required: ');
ak = @(f,x,k) int(f*cos(k*pi*x/l)/l,x,L,U);
bk = @(f,x,k) int(f*sin(k*pi*x/l)/l,x,L,U);
fs = @(f,x,n) ak(f,x,0)/2 + symsum(ak(f,x,k)*cos(k*pi*x/l) + bk(f,x,k)*sin(k*pi*x/l),k,l,n);
fst = ak(f,x,0)/2;
for i=1:n
   fst = fst + ak(f,x,i)*cos(i*pi*x/l) + bk(f,x,i)*sin(i*pi*x/l);
   disp(['Harmonics upto:', num2str(i) ]);
   disp(fst);
  ezplot(f,[L U]);
   hold on;
  ezplot(fst,[L U]);
   %set(h,'LineStyle','-','color',[i/n 1/i 1/n],'LineWidth',2);
   title(['Parallel sums up to n=',num2str(i) ])
end

9.  Legandre

X = linspace(-1,1,500); X = X(:);
N = 10;
Y = zeros(numel(X),N);
Y(:,1) = 1;
Y(:,2) = X;
for n=1:(N-1)
 Y(:,n+2) = ((2*n+1) .* X .* Y(:,n+1) - n .* Y(:,n)) / (n+1);
end
figure('Position',[0 0 1024 768])
19
plot(X, Y(:,1:6), 'LineWidth', 1.5)
axis([-1 1 -1.1 1.1])
grid on;
legend(num2str((1:6)'-1,'P_%d(x)'),'Location','Best')
xlabel('x')
ylabel('P_n(x)')
title('Legendre Polynomials')

10. HARMONIC ANALYSIS

Aim: From the given data points of a function , generate the first few harmonics of Fourier series and visualize them () f x
MATLAB Code:
clear all
clc
syms t
n = input('Enter the number of data points n : ');
x_0 = input('Enter the starting value of x : ');
s = input('Enter the length of the spacing between successive values of x : ');
count = input('Enter 0 if the unit of x is degrees, otherwise enter a non-zero number: ');
y = input('Enter the y values (as a row vector): ');
n1 = input('Enter the number of harmonic of the series n1 : ');
x = x_0 + (0:n-1)*s; if(count == 0)
x=x*pi/180;
s=s*pi/180;
end
l=0.5*(x(n)+s-x(1)); % l=pi if it is degree
a_0 = (2/n)*sum(y);
for i=1:n1
yc=y.*cos(i*pi*x/l);
ys=y.*sin(i*pi*x/l);
a(i)=(2/n)*sum(yc);
b(i)=(2/n)*sum(ys);
end
F_s=a_0/2;
figure('Position',[1 1 1024 768])
for i=1:n1
subplot(n1,1,i)
plot(x,y,'r*');
hold on
F_s = F_s+a(i).*cos(i*pi*t/l)+b(i).*sin(i*pi*t/l);
ezplot(F_s, [x(1) x(n)])
title(['Harmonics up to: ', num2str(i)])
end
disp('Fourier series : ')
vpa(F_s,4)

11. SOLUTION FOR SECOND ORDER DIFFERENCE EQUATIONS WITH CONSTANT COEFFICIENTS

Aim: To write the MATLAB code to find the solution for the second order linear difference equation with constant co-efficient.

MATLAB Code
clear all
clc
syms n k1 k2
a = input('Enter the coefficient of y(n+2): ');
b = input('Enter the coefficient of y(n+1): ');
c = input('Enter the coefficient of y(n): ');
g = input('Enter the non-homogeneous part: ');
r = roots([a b c]);
if imag(r)~=0
rho = sqrt(real(r(1))^2 + imag(r(1))^2);
theta = atan(abs(imag(r(1)))/real(r(1)));
y1 = (rho^n)*cos(n*theta);
y2 = (rho^n)*sin(n*theta);
elseif r(1)==r(2)
y1 = r(1)^n;
y2 = n*r(1)^n;
else
y1 = r(1)^n;
y2 = r(2)^n;
end
Co = det([y1, y2;subs(y1,n,n+1), subs(y2,n,n+1)]); %Casoratian of the solutions
y_c = k1*y1 + k2*y2;
disp('Complementary Solution is: ');
disp(y_c);
if(g ~= 0)
y11 = subs(y1,n,n+1);
y21 = subs(y2,n,n+1);
Co1 = subs(Co,n,n+1);
u1 = simplify(symsum(-g*y21/Co1,n,0,n-1));
u2 = simplify(symsum(g*y11/Co1,n,0,n-1));
y_p = simplify(u1*y1+u2*y2);
y = y_c + y_p;
else
y = y_c;
end
check = input('If the given problem has initial conditions then enter 1 else enter 0: ');
if (check == 1)
yval1 = input('Enter the initial condition at n = 0: ');
yval2 = input('Enter the initial condition at n = 1: ');
cond1 = strcat(char(subs(y,n,0)),'=',num2str(yval1));
cond2 = strcat(char(subs(y,n,1)),'=',num2str(yval2));
[k1,k2] = solve(cond1,cond2);
y = subs(y);
end
disp('Complete Solution is: ')
disp(collect(collect(y,y1),y2))
if(check ~= 0)
nrange = linspace(0,10,40);
Y = subs(y,n,nrange);
stem(nrange,Y);
set(gca,'XTick',linspace(0,10,11))
xlabel('n');
ylabel('y(n)');
end


12. SOLUTION OF FIRST AND SECOND ORDER DIFFERENCE EQUATIONS USING Z - TRANSFORM

Aim: To solve and visualize Difference equation using Z – Transform in MATLAB.
MATLAB CODE for solving second order difference Equations:
clc
clear all
syms n z Yz g
yn = sym('y(n)');
yn1 = sym('y(n+1)');
yn2 = sym('y(n+2)');
a = input('The Coefficient of y(n+2) = ');
b = input('The Coefficient of y(n+1) = ');
c = input('The Coefficient of y(n) = ');
nh = input('The Non Homogeneous part ');
eq = a*yn2 + b*yn1 + c*yn - nh;
Zeq = ztrans(eq, n, z);
if (a==0)
d = input('The initial value at 0 is ');
Zeq = subs(Zeq,{'ztrans(y(n), n, z)', 'y(0)'}, {Yz, d});
else
d = input('The initial value at 0 is ');
e = input('The initial value at 1 is ');
Zeq = subs(Zeq,{'ztrans(y(n), n, z)', 'y(0)', 'y(1)'}, {Yz, d, e});
end
eq = collect(Zeq, Yz);
Y = solve(eq, Yz);
Y = iztrans(Y, z, n);
Y = simple(Y)
m = 0:20;
y = double(subs(Y,n,m));
stem(y)
title('Difference equation');
xlabel('n'); ylabel('y(n)');

13. APPLICATIONS OF SECOND ORDER DIFFERENTIAL EQUATIONS (Springs)

AIM : To know how the second order differential equations can be applied to solve problems concerning the vibrations of springs using MATLAB.

MATLAB Code:
clear
clc
syms x k1 k2
a=input('Enter the coefficient of D2y:');
b=input('Enter the coefficient of Dy:');
c=input('Enter the coefficient of y:');
a1=a/a; b1=b/a; c1=c/a;
eq='a1*x^2+b1*x+c1=0'; %%% solving of homogeneous part
eq1=subs(eq,'[a1,b1,c1]',[a1,b1,c1]);
r=solve(eq1, 'x');
if image(r)~=0
disp('under damping:')
y1=exp(real(r(1))*x)*cos(imag(r(1))*x);
y2=exp(real(r(1))*x)*sin(abs(imag(r(1)))*x);
elseif r(1)==r(2)
disp('critical damping:')
y1=exp(r(1)*x);
y2=x*exp(r(1)*x);
else
disp('over damping:')
y1=exp(r(1)*x);
y2=exp(r(2)*x);
end
W=simplify(y1*diff(y2)-y2*diff(y1)); %%W is the Wronskian
y_c=k1*y1+k2*y2;
disp('The complementary sol y_c is:')
g=input('Enter the non-homogeneous part:'); %%% Solving of Non-homogeneous
y_p=(1/a)*(-y1*int(y2*g/W)+y2*int(y1*g/W));
y=y_c + y_p;
disp('The particular sol y_p is:')
disp(y_p)
disp('hence the general sol y is:')
disp(y)
i_1=input('The initial conditions x and the value of y at x in vector form:')
i_2=input('The initial conditions x and Dy at x in vector form:')
i1=strcat(char(subs(y,x,i_1(1))),'=',num2str(i_1(2)));
i2=strcat(char(subs(diff(y),x,i_2(1))),'=',num2str(i_2(2)));
[k1, k2]=solve(i1,i2);
C_S=vpa(subs(y),3);
disp('Required complete solution is:')
disp(C_S)
ezplot(C_S,[0 1.5])


14. APPLICATIONS OF SECOND ORDER DIFFERENTIAL EQUATIONS (Electric Circuits)

AIM
To know how the second order differential equations can be applied to solve problems concerning the analysis of electric circuits using MATLAB

Same Code as previous

MATLAB Output
In the Command Window:
Enter the coefficient of D2Q: 1
Enter the coefficient of DQ: 40
Enter the coefficient of Q: 625