help with my 4th order runge kutta code

1 回表示 (過去 30 日間)
Vezzaz
Vezzaz 2021 年 10 月 17 日
回答済み: Vezzaz 2021 年 10 月 18 日
I received this error: "Unable to perform assignment because the size of the left side is 1-by-4 and the size of the right side is 4-by-4." when trying to run my 4th order runge kutta code to solve the 4 non linear odes of an elastic pendulum. Any help on how to fix this would be appreciated. The line of code where the error is:
Y(i+1,:) = Y(i,:) + phi*h;
edit: here is my code
tspan1=[0,2*pi];
h=0.01;
r1=11;
rdot1=0;
theta1=0;
thetadot1=0;
init1= [r1;theta1;rdot1;thetadot1];
[T,Y]=rk4_fun(@spring1,tspan1,init1,h);
function [T,Y]=rk4_fun(f,tspan,y0,h)
x0 = tspan(1);xf = tspan(2);
T = (x0:h:xf)';
n = length(T);
n_eqn = length(y0);
Y = ones(n,n_eqn);
Y(1,:) = y0;
for i = 1:n-1
k1 = f(T(i),Y(i,:))';
k2 = f(T(i)+h/2,Y(i,:) + k1*h/2)';
k3 = f(T(i)+h/2,Y(i,:) + k2*h/2)';
k4 = f(T(i)+h,Y(i,:) + k3*h)';
phi = (k1+2*k2+2*k3+k4)/6;
Y(i+1,:) = Y(i,:) + phi*h;
end
end
function ydt = spring1(t, y)
g = 1000;
km = 100;
lo=10;
[rows,cols] = size(y);
ydt = zeros(rows,cols);
ydt(1) = y(3);
ydt(2) = y(4);
ydt(3) = y(1)*y(4)*y(4) + g*cos(y(2)) - km*(y(1)-lo);
ydt(4) = -(g*sin(y(2))+2*y(3)*y(4))/y(1);
end
  2 件のコメント
James Tursa
James Tursa 2021 年 10 月 17 日
In the future, it is best to post the entire code that produces the error so that we can get the proper context. Without seeing the rest of your code we have to make guesses.
John D'Errico
John D'Errico 2021 年 10 月 17 日
Without seeing the rest of the code, it is impossble to help you.
Most importantly, we need to know the size of both the phi and h variables. Not what you THINK them to be. But what they are. EXACTLY.

サインインしてコメントする。

採用された回答

James Tursa
James Tursa 2021 年 10 月 17 日
My guess without seeing the rest of your code is that phi is a column vector, so when you add phi*h to the row vector Y(i,:) it creates a 4x4 matrix. Try changing phi to a row vector, or transpose it in the equation itself. E.g.,
Y(i+1,:) = Y(i,:) + phi'*h;

その他の回答 (1 件)

Vezzaz
Vezzaz 2021 年 10 月 18 日
i posted the entire code for this part.

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by