Runge-Kutta 4th order function error (Matrix dimensions must agree)
1 view (last 30 days)
Show older comments
hello, i'm trying to make a Runge-Kutta 4th order function, and it works for some questions but i can't get it to work for this one, please check my code
close all; clc; clear all;
tspan = [0 10];
xi = [1;1];
params = [1 1]; %a = 1, b = 1 ; see function file
[t,x] = ode45(@hw3p4_func,tspan,xi,[],params);
y0 = x(1); y1 = x(2); ti = tspan(1) ; tf = tspan(2);
[t_rk, x_rk] = RK_4(@hw3p4_func, xi,y0,0.1);
plot(t,x(:,1))
%function to be solved
function dxdt = hw3p4_func(t,x, params)
y1 = x(1);
y2 = x(2);
a = params(1);
b = params(2);
dxdt = [a*y2;
(-b*y1+(1-y1^2)*y2)];
end
%RK code
function [x,y] = RK_4(func, x,y0,delta_x)
y = zeros(1,length(x));
y(1) = y0;
n = length(x)-1;
for i = 1:n
k1 = func(x(i),y(i));
k2 = func(x(i)+.5.*delta_x,y(i)+.5.*k1.*delta_x);
k3 = func(x(i)+.5.*delta_x,y(i)+.5.*k2.*delta_x);
k4 = func(x(i)+ delta_x,y(i)+ k3.*delta_x);
y(i+1) = y(i)+((k1+2.*k2+2.*k3+k4)./6).*delta_x;
end
0 Comments
Answers (1)
James Tursa
on 8 Jul 2020
Edited: James Tursa
on 8 Jul 2020
Your RK_4 function is not set up to handle vector equations ... it is only set up to handle scalar equations. Also you are not passing params into RK_4, so the derivative function will not have that needed input.
So when you call it, you need to pass in an initial column vector and params in the function handle. E.g.,
[t_rk, x_rk] = RK_4(@(t,x)hw3p4_func(t,x,params), tspan, xi,0.1);
And the derivative function fixes would look something like this:
function [x,y] = RK_4(func, tspan,y0,delta_x)
n = round((tspan(2)-tspan(1))/delta_x)+1; % new calculation of n
x = linspace(tspan(1),tspan(2),n); % create a vector of x's from the endpoints with n elements
y = zeros(numel(y0),n); % changed to matrix of column vectors
y(:,1) = y0; % changed y to column vector
for i = 1:n-1 % changed top index on rhs to n-1
k1 = func(x(i) ,y(:,i) ); % changed y to column vector
k2 = func(x(i)+.5.*delta_x,y(:,i)+.5.*k1.*delta_x); % changed y to column vector
k3 = func(x(i)+.5.*delta_x,y(:,i)+.5.*k2.*delta_x); % changed y to column vector
k4 = func(x(i)+ delta_x,y(:,i)+ k3.*delta_x); % changed y to column vector
y(:,i+1) = y(:,i)+((k1+2.*k2+2.*k3+k4)./6).*delta_x; % changed y to column vector
end
end
You could plot the results as (picking off 1st row since states are stored as column vectors):
plot(t_rk,x_rk(1,:));
0 Comments
See Also
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!