Runge-Kutta 4th order function error (Matrix dimensions must agree)

2 visualizaciones (últimos 30 días)
ragheed idrees
ragheed idrees el 8 de Jul. de 2020
Editada: James Tursa el 8 de Jul. de 2020
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

Respuestas (1)

James Tursa
James Tursa el 8 de Jul. de 2020
Editada: James Tursa el 8 de Jul. de 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,:));

Categorías

Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.

Productos


Versión

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by