Project 1 - Data Handling and Interpolation

Project 1:

Spring - Mass System

% mass = 10 kg

m = input('Enter the value of mass (in kgs): ');

% spring stiffness = 11 N/mm

k = input('Enter the spring constant (in N/mm): ');

t_span = [0 15]; % Time span

x_ic = [8 0]; % IC: x(1) for displacement and x(2) for velocity

spring_mass_func = @(t,x) [x(2); -(k/m)*x(1)];

% first row integrates velocity and 2nd row integrates acceleration

%since its is a numerical integration, the output from the integration wil

% again act as an input.

% Solve using ODE23

[tout1,xout1] = ode23(@(t,x) spring_mass_func(t,x),t_span,x_ic);

x_dot1 = xout1(:,2); %velocity values

x1 = xout1(:,1); % displacement values

%% solve for acceleration

x_dot_dot1 = -(k/m)*x1;

figure(1)

plot(tout1,x1,'b--'); %displacement plot

hold on;

plot(tout1,x_dot1,'g--',tout1,x_dot_dot1,'r-'); %velcoity & acceleration plot

xlabel('Time (in sec)');

ylabel({'x (in mm), $\dot{X}$ (in $\mathrm{mm/s}$), $\ddot{X}$ (in $\mathrm{mm/s}^2$)'},'Interpreter','latex');

title('Solution of ODE23')

legend({'X','$\dot{X}$','$\ddot{X}$'},'Interpreter', 'latex');

hold off

% Solve using ODE45

[tout2,xout2] = ode45(@(t,x) spring_mass_func(t,x),t_span,x_ic);

x_dot2 = xout2(:,2); %velocity values

x2 = xout2(:,1); % displacement values

%% solve for acceleration

x_dot_dot2 = -(k/m)*x2;

figure(2)

plot(tout2,x2,'c--'); %displacement plot

hold on;

plot(tout2,x_dot2,'m--',tout2,x_dot_dot2,'k-'); %velcoity & acceleration plot

xlabel('Time (in sec)');

ylabel({'x (in mm), $\dot{X}$ (in $\mathrm{mm/s}$), $\ddot{X}$ (in $\mathrm{mm/s}^2$)'},'Interpreter','latex');

title('Solution of ODE45')

legend({'X','$\dot{X}$','$\ddot{X}$'},'Interpreter', 'latex');

hold off

%compare the solutions of ODE solvers

figure(3)

% Displacement vs. Time

subplot(3, 1, 1);

plot(tout1, x1, 'r-', tout2, x2, 'm-');

xlabel('Time (in sec)');

ylabel('x (in mm)');

legend('ODE23', 'ODE45');

title('Displacement vs. Time');

% Velocity vs. Time

subplot(3, 1, 2);

plot(tout1, x_dot1, 'r-', tout2, x_dot2, 'm-');

xlabel('Time (in sec)');

ylabel({'$\dot{x}$ (in $\mathrm{mm/s}$)'},'Interpreter', 'latex');

legend('ODE23', 'ODE45');

title('Velocity vs. Time');

% Acceleration vs. Time

subplot(3, 1, 3);

plot(tout1, x_dot_dot1, 'r-', tout2, x_dot_dot2, 'm-');

xlabel('Time (in sec)');

ylabel({'$\dot{x}$ (in $\mathrm{mm/s}^2$)'}, 'Interpreter', 'latex');

legend('ODE23', 'ODE45');

title('Acceleration vs. Time');

Observation:

Usually problems of this type are not difficult to solve by both ODE23 and ODE4.

By time stamping I understood the elapsed time for calculation by both solvers took same time.

ODE45 takes longer time to solve usually because it is uses runge-kutta 4th and 5th order for solving differential equation.

Project 3:

Fmin and Max search

%Develop the find function

func = @(x) sin(x.^2)+cos(x) - 2*sin(x).* cos(x);

% Generator Plot

x = linspace(-3*pi,3*pi,200);

y = func(x);

plot(x,y,'r-');

title('Fmin and Global max search');

xlabel('x');ylabel('y');

legend('y');

print('y','-dpng')

%find all the local minima in the range of -3*pi to 3*pi

% since fminsearch func always gives the nearest minima, we have to be sure

% to avoid the repeated minima outputs for given guess input x.

local_minima_x = [];

for i = 1:length(x)-1

k1 = x(i);

k2 = x(i+1);

minima_x(1) = fminsearch(@(x) func(x),k1);

minima_x(2) = fminsearch(@(x) func(x),k1);

if minima_x(1) >= k1 && minima_x(1) <= k2

local_minima_x = [local_minima_x minima_x(1)];

end

end

local_minima_y = func(local_minima_x);

%find the global maximum

[y_val,y_index] = max(y);

global_max_x = x(y_index);

global_max_y = y(y_index);

%High lighting local min and max through text

hold on

plot(local_minima_x,local_minima_y,'go');

for i = 1:length(local_minima_y)

text(local_minima_x(i)+0.05,local_minima_y(i)-0.1, sprintf('Minima: x = %.2f, y = %.2f', local_minima_x(i), local_minima_y(i)), 'FontSize', 5);

end

plot(global_max_x,global_max_y,'c*');

text(global_max_x,global_max_y, sprintf('G_max: x = %.2f, y = %.2f', global_max_x, global_max_y), 'FontSize', 5);

hold off

legend('y','local minima','global maxima');

print('y with minima and maxima','-dpng');

Project 4:

contour plots are similar to mesh plots, takes two independent variables (x and y) and graphs the output (z).

while the mesh plot shows the graph with z axis for all its range (3D plot), the contour represents the same in 2D with z values in isolines, i.e. the intensity of the color of the line represents the z-axis values. (refer the color bar in attached figure).

For example I have simulated z = y*exp(-(x^2+y^2)).

Matlab commands:

Contour plot example

%independent variables define

X = -pi:pi/10:pi;

Y = -pi:pi/10:pi;

%convert to meshgrid for data points

[x,y] = meshgrid(X,Y);

% dependent variable function definition

z = y.*exp(-(x.^2+y.^2));

figure(1)

contour(x,y,z,'ShowText','on');

colorbar;

title('contour plot of z = y*exp(-(x^2+y^2))');

xlabel('x');ylabel('y');

print('contour_ex','-dpng');

figure(2)

mesh(x,y,z);

title('mesh plot of z = y*exp(-(x^2+y^2))');

xlabel('x');ylabel('y');zlabel('z')

print('mesh_ex','-dpng');