Question

2nd order ODE of modeling a cylinder oscillating in still water wate wate Figure 1. A cylinder oscillating in still water. A cylinder floating in the water can be modeled by the second order ODE: dy dy dt dt where y is the distance from the water level to the neutral line as shown in Figure 1. The mass, damping constant, and buoyancy constant of this system are given as: m lkg, c-0.1kg/sec, k-0.01kg/s2. The cylinder rested in the water and suddenly was risen and dropped from y-0.2m. Hence, the initial conditions are y(0)-0.2 m and y'(0)-0 m/sec. /p (0) Use the knowledge you learned in the lecture to develop your Matlab code for Explicit-Euler's method, midpoint method, and the 4th order Runge-Kutta method. Use these methods to predict the y() from t-0(sec) to t-200(sec) and compare your results with the exact solution of c-0.1kg/s, m-l kg, k-0.01kg/s2. which is given as: y(1) = exp(-0.051)[02cos(at) +0.1 1547sin(at)] where ω = 0.086602s-ad/s Please use the time step of 0.1(sec), 0.5(sec), 1(sec), 5(sec), 10(sec) to run your numerical solutions and discuss your findings.

I mostly needed help with developing matlab code using the Euler method to create a graph. All the other methods are doable once I have a proper Euler method code to refer to.

Answer 1

The Euler Method in Matlab For the Second order DE is given below.There are two files one is function file which is contains the euler method. And the another is script file that prepares the differential equation and sets the constants and then feeds the DE into euler function

************************************euler_method.m**********************

function [t,y] = euler_method( deriv_func,t0,dt,t_end,y0 )
%deriv_func is the differential eq.
%t0, t_end are the time range. dt is the increment
%calculate the size of output array
size=ceil((t_end-t0)/dt);

time_range=t0:dt:t_end;

%following line creates array with all elements 0
Y=zeros(size+1,2);
%initialize the Y with initial values
Y(1,:)=y0;

%the below loop solves the differential equation
% y(i+1)=y(i)+h*f(t,y)
for i=1:size
%x will give us column vector with two values
%one for y(1) and another for y(2)
x=deriv_func(time_range(i),Y(i,:));
%below is the euler method formula
Y(i+1,:)=Y(i,:)+dt*x';

end
t=time_range;
y=Y;
end

****************************SecondOrderDE.m************************

clc
clear all;

m=1;
c=0.1;
k=0.01;

initial=[0.2,0];
%Here we convert the Second order
%ODE in to first order by renaming the
%variables. y'=y(2).So differential of y' i.e. y"
% is d(y(2))/dt
%Now d(y(2))/dt = (-k*y(1)-c*y(2))/m.
%That is what we are representing below using Anonymous function
F= @(t,y)[y(2);(-k*y(1)-c*y(2))/m];

t0=0;
t_end=100;
dt=0.001;
[t,y]=euler_method(F,t0,dt,t_end,initial);

%Following is the Analytical function against which we compare
%Euler method.
Y=@(t)exp(-0.05.*t).*(0.2.*cos(0.0866025.*t)+0.11547.*sin(0.0866025.*t))
y_analytical=Y(t);

%Now we compare euler with True method
for i=1:length(t)
fprintf('Euler:%f\t\tAnalytical:%f\n',y(i,1),y_analytical(i))
end

%plot (t,y(;,1)) plots displacement vs time

plot(t,y(:,1)), xlabel('time'), ylabel('Displacement'),title('distance-time');

Euler:-0.000433 Euler:-0.000433 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000434 Euler:-0.000435 Euler:-0.000435 Euler:-0.000435 Analytical:-0.000433 Analytical:-0.00043:3 Analytical:-0.00043:3 Analytical:-0.00043:3 Analytical:-0.000433 Analytical:-0.00043:3 Analytical:-0.00043:3 Analytical:-0.00043:3 Analytical:-0.000433 Analytical:-0.000434 Analytical:-0.000434 Analytical:-0.000434 Analytical:-0.000434 Analytical:-0.000434

alstance-time 0.2 0.15 50.1 CO 0.05 0 0.05 0 10 20 304050 60 70 80 90 100 time