Runge-Kutta 4th order function error (Matrix dimensions must agree)
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
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,:));
0 Kommentare
Siehe auch
Kategorien
Mehr zu Get Started with MATLAB finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!