% March 17
%% Example 1, use ode45 to solve for a system of 3 ODEs
clear,clc
close all
[t,y] = ode45(@dydtfun,[0,12],[0,1,1]);
plot(t,y(:,1),t,y(:,2),t,y(:,3))
xlabel('time')
ylabel('y values')
legend('y1','y2','y3')
%% Example 2, use ode45 to solve for a vibration system
y0 = [0.5,0];
tspan = [0,50];
[t,y] = ode45(@vibfun,tspan,y0);
figure
subplot(2,1,1)
plot(t,y(:,1))
subplot(2,1,2)
plot(t,y(:,2))
%% Example 2 modified, compare ode45 and ode15s
clear,clc
close all
y0 = [0.5, 0];
tspan = [0,1e6]; % 1e6 is 10^6
tic
[t,y] = ode45(@vibsfun,tspan,y0);
toc
tic
[t1,y1] = ode15s(@vibsfun,tspan,y0);
toc
plot(t,y(:,1),t,y(:,2),t1,y1(:,1),'o',t1,y1(:,2),'o')
legend('x by ode45','v by ode45','x by ode15s','v by ode15s')
%% Define all the ode functions
% Example 1
function yprime = dydtfun(t,y)
yprime = [0;0;0];
yprime(1) = y(2)*y(3)*t;
yprime(2) = -y(1)*y(3);
yprime(3) = -0.51*y(1)*y(2);
end
% alternatively
% yprime = [y(2)*y(3)*t;-y(1)*y(3);-0.51*y(1)*y(2)]
% Example 2
function yprime = vibfun(t,y)
% y(1) is x
% y(2) is xdot
% yprime(1) is xdot or y(2)
% yprime(2) is x2dot
m = 10;
k = 4;
c = 2;
f0 = 0.05;
w = 2;
yprime = [0;0];
yprime(1) = y(2);
yprime(2) = f0/m*sin(w*t) - c/m*y(2) - k/m*y(1);
end
% Example 2 modified to be a stiff ode problem
function yprime = vibsfun(t,y)
% y(1) is x
% y(2) is xdot
% yprime(1) is xdot or y(2)
% yprime(2) is x2dot
m = 1;
k = 1e-3;
c = 1;
yprime = [0;0];
yprime(1) = y(2);
yprime(2) = - c/m*y(2) - k/m*y(1);
end
% Mar. 22
%% Beam deflection
% define all the known parameters
% divide the length into N subdivisions
%
N = 10000;
dx = L/N;
N1 = round(L1/dx);
N2 = round((L2-L1)/dx);
% calculate the reactions forces
R1 =
R2 =
% create the known vector, p_n
Pn = zeros(N+1,1);
% replace the elements in Pn using the specific expression/values for that
% beam section
for n = 2:N1+1
Pn(n) =
end
for n = N1+2:N1+N2+1
Pn(n) =
end
for n = N1+N2+2:N
Pn(n) =
end
% define the coefficient matrix A
A = zeros(N+1,N+1);
% use diag function to modify the matrix A, matching the tridiagnol matrix
% solve for the linear equations
y = A\Pn;
