Info

This question is closed. Reopen it to edit or answer.

How can i implement the damping coefficient into my codes?

1 view (last 30 days)
clear all
clc
c = [5 40 200];
m = 20;
k = 20;
h = 0.5;
tn = 15;
y = [1 0];
y01 = y(1);
y02 = y(2);
[t,x,y] = RK2516(@dy1,@dy2,tn,y01,y02,h,c,m,k);
disp([t,x,y]);
----------------------
function dy = dy1(t,y1,y2,c,m,k)
dy = y2;
--------------
function dy = dy2(t,y1,y2,c,m,k)
dy = -(c*y2 - k*y1)/m;
--------------------------
function [t,y1,y2] = RK2516(dy1,dy2,tn,y01,y02,h,c,m,k)
% dy1 = m-file that evaluates the first ODE
% dy2 = m-file that evaluates the second ODE
% tn = range of initial t value to final t value
% y01 = initial value of the first dependent variable
% y02 = initial value of the second dependent variable
% h = step size
t = [0:h:tn]';
n = length(t);
h2 = h/2;
h3 = h/3;
h6 = h/6;
y1 = y01*ones(n,1);
y2 = y02*ones(n,1);
if t(n) < tn
t(n+1) = tn;
n = n + 1;
end
for j = 2:n
k1 = feval(dy1, t(j-1), y1(j-1), y2(j-1));
q1 = feval(dy2, t(j-1), y1(j-1), y2(j-1));
k2 = feval(dy1, t(j-1) + h2, y1(j-1) + h2*k1, y2(j-1) + h2*q1);
q2 = feval(dy2, t(j-1) + h2, y1(j-1) + h2*k1, y2(j-1) + h2*q1);
k3 = feval(dy1, t(j-1) + h2, y1(j-1) + h2*k2, y2(j-1) + h2*q2);
q3 = feval(dy2, t(j-1) + h2, y1(j-1) + h2*k2, y2(j-1) + h2*q2);
k4 = feval(dy1, t(j-1) + h, y1(j-1) + h*k3, y2(j-1) + h*q3);
q4 = feval(dy2, t(j-1) + h, y1(j-1) + h*k3, y2(j-1) + h*q3);
y1(j) = y1(j-1) + h6*(k1 + k4) + h3*(k2 + k3);
y2(j) = y2(j-1) + h6*(q1 + q4) + h3*(q2 + q3);
end
end

Answers (0)

Community Treasure Hunt

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

Start Hunting!