Using Runge-Kutta algorithm to solve second order ODE

27 次查看（过去 30 天）
Daniel Jordan
Daniel Jordan 2022-4-22
编辑： Torsten 2022-4-22
For this task I need to create a program which manually uses the RK algorithm to solve a second order ODE. I'm trying to get my head around it but I think I'm going at it in the wrong way.
My original function is with given values for m, F0 and w, so I have rearranged it to x'' = 50sin(8*pi*t)-x'-50x
I've then switched x' for y and rearranged to get 50sin(8*pi*t) = y' + y + 50x and this is what I'm using in the script. I need to plot x(t) and y(t) but at this point I'm just trying to find values for x and y over time.
The method I've used to get this far is below the code. Excuse my poor code. Thanks in advance.
PS I'm not allowed to use a function like ode45 to solve it for me.
close all; clear; clc;
t=0;x=0;y=0;
dy=myfunction(t,x,y)
x0=0;y0=0
h = 0.5;
for t=0:5.5
dx1 = h*y
dy1 = h*dy
dx2 = h*(y+(dy1/2))
dy2 = h*myfunction(t+(h/2),x+(dx1/2),y+(dy1/2))
dx3 = h*(y+(dy2/2))
dy3 = h*myfunction(t+(h/2),x+(dx2/2),y+(dy1/2))
dx4 = h*(y+dy3)
dy4 = h*myfunction(t+h,x+dx3,y+dy1)
deltax = (dx1+(2*dx2)+(2*dx3)+dx4)/6
deltay = (dy1+(2*dy2)+(2*dy3)+dy4)/6
xth = x+deltax
yth = y+deltay
x=xth
y=yth
t = t+h
end
function [dy] = myfunction(t,x,y)
dy=(50*sin(8*pi*t))-(50*x)-y;
end
Method used (I swapped v for y):

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2022-4-22
Can you take it from here ?
%Solves y'' - (-exp(-B*t)-y+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t) = 0
% t in [0 1]
% y(0) = 1, y'(0) = -1
% B = 4
tstart = 0.0;
tend = 1.0;
h = (tend - tstart)/20;
T = (tstart:h:tend).';
Y0 = [1 -1];
B = 4;
f = @(t,y) [y(2) -exp(-B*t)-y(1)+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t];
Y = runge_kutta_RK4(f,T,Y0);
plot(T,Y)
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
  2 个评论
Daniel Jordan
Daniel Jordan 2022-4-22
So can I change this line
f = @(t,y) [y(2) -exp(-B*t)-y(1)+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t];
into my own equation? Wouldn't it need to have t, x and y?
Also what do y(1) and y(2) mean in that line?
Torsten
Torsten 2022-4-22
编辑：Torsten 2022-4-22
For your function,
f = @(t,y) [y(2) F0/m*sin(omega*t)-c/m*y(2)-k/m*y(1)]
where
y(1) = y, y(2) = y'.
So you shouldn't work with x and y, but with a vector y=(y(1),y(2)) to make everything easier.
Of course, omega, F0, m and k have to be given values before.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by