Project 1 - Data Handling and Interpolation
Project 1
% Simple Harmonic Oscillator Solver
% Get mass (m) and spring constant (k) from the user
prompt = {'Enter mass (m):', 'Enter spring constant (k):'};
dlgtitle = 'Input';
dims = [1 35];
definput = {'1', '1'};
user_input = inputdlg(prompt, dlgtitle, dims, definput);
% Convert user input to numeric values
m = str2double(user_input{1});
k = str2double(user_input{2});
% Define the differential equation for simple harmonic oscillator
ode_eq = @(t, y) [y(2); -k/m * y(1)];
% Initial conditions: x(0) = 1, x_dot(0) = 0
initial_conditions = [1; 0];
% Solve using ode45
[t_ode45, sol_ode45] = ode45(ode_eq, [0 10], initial_conditions);
% Solve using ode23
[t_ode23, sol_ode23] = ode23(ode_eq, [0 10], initial_conditions);
% Plot results
figure;
% Plot displacement (x)
subplot(3,1,1);
plot(t_ode45, sol_ode45(:,1), '-o', t_ode23, sol_ode23(:,1), '-s');
title('Displacement (x)');
legend('ode45', 'ode23');
% Plot velocity (x_dot)
subplot(3,1,2);
plot(t_ode45, sol_ode45(:,2), '-o', t_ode23, sol_ode23(:,2), '-s');
title('Velocity (ẋ_dot)');
legend('ode45', 'ode23');
% Plot acceleration (x_dot_dot)
subplot(3,1,3);
acc_ode45 = -k/m * sol_ode45(:,1);
acc_ode23 = -k/m * sol_ode23(:,1);
plot(t_ode45, acc_ode45, '-o', t_ode23, acc_ode23, '-s');
title('Acceleration (ẍ_dot_dot)');
legend('ode45', 'ode23');
% Adjust plot layout
xlabel('Time');
sgtitle('Simple Harmonic Oscillator');
% Explanation
disp('Comparison:');
disp('ode45 generally provides a more accurate solution for stiff ODEs,');
disp('while ode23 is often more efficient for non-stiff problems.');
disp('In this case, since the harmonic oscillator is not stiff,');
disp('ode45 and ode23 should produce similar results.');
simple harmonic oscillator equation using the 'ode45' and 'ode23' solvers. It also plots the displacement (x), velocity (ẋ), and acceleration (ẍ) over time. The code uses the
inputdlg function to get 'm' and 'k' values from the user.In this code, the user is prompted to enter values for mass (
m) and spring constant (k). The simple harmonic oscillator equation is then solved using both ode45 and ode23 solvers. The results are plotted for displacement (x), velocity (ẋ), and acceleration (ẍ). The explanation at the end clarifies that for non-stiff problems like the harmonic oscillator, ode45 and ode23 should produce similar results.Plot Comparison:
- When you compare the plots for displacement (x), velocity (ẋ), and acceleration (ẍ) obtained from
ode45andode23, you'll likely observe that they closely overlap. - The differences, if any, might be more apparent in the error behavior or the number of function evaluations performed by each solver.
ode45is often more accurate thanode23for smooth solutions and non-stiff problems. It uses a variable step size, adaptive Runge-Kutta method (specifically the Dormand-Prince method) that adjusts the step size to control the local error.ode23uses a lower-order, fixed step size method (Bogacki-Shampine method). While it is efficient for non-stiff problems, it may require smaller step sizes to achieve the same level of accuracy asode45.
ode23is typically more efficient for non-stiff problems because it uses a lower-order method with fewer function evaluations per step. It may converge faster when the solution is relatively smooth.ode45may be more computationally expensive due to its higher-order method and variable step size, especially if the solution exhibits rapid changes.
Project 2
[data1] = readmatrix('isa.xlsx'); % loading given Temperature and Altitude excel data sheet using readmatrix command.
H = data1(:, 1); % consider H is altitude.
T = data1(:, 2); % consider T is a Temperature.
p = polyfit(H, T, 3); % use of polyfit function to determine regression model T=f(h).
figure(1);
plot(H, T);
hold on;
xlabel('Altitude (m)'); % labeling the figure
ylabel('Temperature (\circC)');
title('ISA plot'); % give title to the figure.
y1 = polyval(p, H); % evaluate the model using polyval function.
plot(H, y1, 'r--');
% Store polyval evaluated data into a text file
evaluatedData = [H, y1];
dlmwrite('evaluated_data.txt', evaluatedData, 'delimiter', '\t', 'precision', 6);
hold off; % Release the hold on the plot
% Additional code if you want to display the coefficients of the polynomial
fprintf('Polynomial Coefficients: %.6f %.6f %.6f %.6f\n', p);
% Additional code if you want to display the equation of the polynomial
equation = sprintf('T = %.6fH^3 + %.6fH^2 + %.6fH + %.6f', p(1), p(2), p(3), p(4));
disp(equation);
- store polyval evaluated data into a text file.
2000 307.536
-1500 302.543
-1000 297.711
-500 293.039
0 288.525
500 284.167
1000 279.962
1500 275.91
2000 272.007
2500 268.253
3000 264.644
3500 261.179
4000 257.855
4500 254.672
5000 251.626
5500 248.716
6000 245.94
6500 243.296
7000 240.781
7500 238.394
8000 236.133
8500 233.995
9000 231.979
9500 230.083
10000 228.305
10500 226.642
11000 225.093
11500 223.655
12000 222.327
12500 221.107
13000 219.992
13500 218.981
14000 218.072
14500 217.262
15000 216.549
15500 215.933
16000 215.41
16500 214.978
17000 214.636
17500 214.381
18000 214.212
18500 214.126
19000 214.122
19500 214.197
20000 214.35
22000 215.696
24000 218.12
26000 221.494
28000 225.689
30000 230.577
Project 3
% Define the function
f = @(x) sin(x.^2) + cos(x) - 2*sin(x).*cos(x);
% Generate a range of x values for plotting
x_values = linspace(-2*pi, 2*pi, 1000);
% Evaluate the function at each x value
y_values = f(x_values);
% Plot the function
figure;
plot(x_values, y_values, 'LineWidth', 1.5);
title('Plot of y = sin(x^2) + cos(x) - 2sin(x)cos(x)');
xlabel('x');
ylabel('y');
% Find local minimum using fminsearch
local_minima_x = fminsearch(f, 0);
local_minima_y = f(local_minima_x);
% Find global maximum using max command
[global_maxima_y, global_maxima_index] = max(y_values);
global_maxima_x = x_values(global_maxima_index);
% Highlight local minima and global maxima on the plot using gtext
hold on;
gtext(['Local Minima: x = ', num2str(local_minima_x), ', y = ', num2str(local_minima_y)]);
gtext(['Global Maxima: x = ', num2str(global_maxima_x), ', y = ', num2str(global_maxima_y)]);
hold off;
Project 4
Give a brief note on the “contour” command. Illustrate with an example.
A contour plot is a graphical method to visualize the 3-D surface by plotting constant Z slices called contours in a 2-D format. The contour plot is an alternative to a 3-D surface plot
The contour plot is formed by:
- Vertical axis: Independent variable 2
- Horizontal axis: Independent variable 1
- Lines: iso-response values, can be calculated with the help (x,y).
The independent variable usually restricted to a regular grid. The actual techniques for determining the correct iso-response values are rather complex and almost always computer-generated.
The contour plot is used to depict the change in Z values as compared to X and Y values. If the data (or function) do not form a regular grid, you typically need to perform a 2-D interpolation to form a regular grid.
For one variable data, a run sequence/ histogram is considered necessary. For two-variable data, a scatter plot is considered necessary. The contour plots can also polar co-ordinates (r,theta) instead of traditional rectangular (x, y, z) coordinates.
Types of Contour Plot:
- Rectangular Contour plot: A projection of 2D-plot in 2D-rectangular canvas. It is the most common form of the contour plot.
- Polar contour plot: Polar contour plot is plotted by using the polar coordinates r and theta. The response variable here is the collection of values generated while passing r and theta into the given function, where r is the distance from origin and theta is the angle from the positive x axis.
- Ternary contour plot: Ternary contour plot is used to represent the relationship between 3 explanatory variables and the response variable in the form of a filled triangle.
Contour plot can be plotted in different programming languages:
- Python/ Matplotlib: Contour plot can be plotted using plt.contour or plt.contourf functions, where plt is matplotlib.pyplot. The difference between these two that plot.contour generates hollow contour plot, the plt.contourf generated filled.
- Matlab: functions such as contourf (2d-plot) and contour3 (3D-contour) can be used for contour plotting
- R: can create a contour plot with filled.contour functions in R.
x = -2:0.2:2;
y = -2:0.2:3;
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
contour(X,Y,Z,'ShowText','on')