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

2 次查看(过去 30 天)
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

回答(1 个)

James Tursa
James Tursa 2020-7-8
编辑:James Tursa 2020-7-8
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,:));

类别

Help CenterFile Exchange 中查找有关 Get Started with MATLAB 的更多信息

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by