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

2 Ansichten (letzte 30 Tage)
ragheed idrees
ragheed idrees am 8 Jul. 2020
Bearbeitet: James Tursa am 8 Jul. 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

Antworten (1)

James Tursa
James Tursa am 8 Jul. 2020
Bearbeitet: James Tursa am 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,:));

Kategorien

Mehr zu Get Started with MATLAB finden Sie in Help Center und File Exchange

Produkte


Version

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by