Using MatLAB to solve 2nd Order DDE, weird output graph but not sure why?

3 views (last 30 days)
Hi, I have an assignment for Uni to solve the 2nd Order ODE:
m*x''(t) + c*x'(t) + kx = P(t)
given initial conditions, x(t=0) = 1.0, dx/dt(t=0) = 0.0, and the values of parameters as m=1 kg, c=0.1 N•s/m, k=42 N/m , P(t)=0 N, and between 0: t : 10, and dt=0.2s,
I have a code which gives me a graph that is incorrect, however I am very new to using MATLAB and I am probably missing something obvious in my code, so any help would be very appreciated.
Code is attatched below!
Thanks!
clc; %Clear Windows%
clear all;
f_1= @(t,V,x) -0.1*V-42*x; %Define Functions%
f_2= @(t,V,x) V;
h=0.2; %Define Time Step and Time Range%
t=0:h:10;
x = zeros(1,length(t)); %Create Matricies for loop%
V = zeros(1,length(t));
x(1)=1; %Input initial conditions%
V(1)=0;
for i=1:50 %Computational loop%
k_1_1 = h*f_1( t(i), x(i), V(i));
k_1_2 = h*f_2( t(i), x(i), V(i));
k_2_1 = h*f_1( t(i)+h/2 , x(i)+k_1_1/2 , V(i)+k_1_2/2 );
k_2_2 = h*f_2( t(i)+h/2 , x(i)+k_1_1/2 , V(i)+k_1_2/2 );
k_3_1 = h*f_1( t(i)+h/2 , x(i)+k_2_1/2 , V(i)+k_2_2/2 );
k_3_2 = h*f_2( t(i)+h/2 , x(i)+k_2_1/2 , V(i)+k_2_2/2 );
k_4_1 = h*f_1( t(i)+h , x(i)+k_3_1 , V(i)+k_3_2 );
k_4_2 = h*f_2( t(i)+h , x(i)+k_3_1 , V(i)+k_3_2 );
xn(i+1)= x(i) + (1/6)*(k_1_1+2*k_2_1+2*k_3_1+k_4_1)*h;
x(i+1) = xn(i+1)+ (1/6)*(k_1_1+2*k_2_1+2*k_3_1+k_4_1)*h;
end
%Exact Solution%
%%%%%%%%%%%%%%%%%%%%%
A=1;
B=0.0077154;
t = 0:0.2:10;
y=A.*exp(-0.05.*t) .* cos(6.48055.*t)+B.*exp(-0.05.*t) .* sin(6.48055.*t);
%%%%%%%%%%%%%%%%%%%%%
%Graphs%
figure1 = figure;
% Create axes
axes1 = axes('Parent',figure1,'FontSize',20,'FontName','Times New Roman');
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');
axis([axes1],[0 1 0 11]);
% Create xlabel
xlabel('t','FontSize',20,'FontName','Times New Roman');
% Create ylabel
ylabel('x','FontSize',20,'FontName','Times New Roman');
% Create multiple lines using matrix input to plot
plot1 = plot(x,t,'Parent',axes1,'LineWidth',2);
plot2 = plot(y,t,'Parent',axes1,'LineWidth',2);
%Label%
set(plot1(1),'Marker','diamond','DisplayName','RK4');
set(plot2(1),'Marker','o','Color',[1 0 0],'DisplayName','Exact');
% Create legend
legend(axes1,'show');

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!