Unable to perform assignment because the size of the left side is 1-by-2 and the size of the right side is 2-by-2. (Line 32)

6 ビュー (過去 30 日間)
function [ t,y ] = rk4sys( f,tspan,y0,h )
%function [ t,y ] = rk4sys( f,tspan,y0,h )
%Uses RK4 method to solve the system ode y1'=p(t,y1,y2) y1(t0)=y1_0
%y2'=q(t,y1,y2) y2(t0)=y2_0
%on the interval tspan=[t0, tend] using time-step h.
%f=[p;q] and is formatted as required by ode45
%y0=[y1_0, y2_0]
%y is the output vector approximating [y1(t), y2(t)] at times stored in t.
%calculate n the number of time-steps needed to reach tend.
n=ceil((tspan(2)-tspan(1))/h);
%Store initial conditions in t and y.
y(1,1)=y0(1);
y(1,2)=y0(2);
t(1)=tspan(1);
%RK4 Iteration Loop
for i=1:n
%If on last time-step, redefine h to make final t value = tend.
if i==n
h=tspan(2)-t(i);
end
k1=f(t(i),y(i,:))';
k2=f(t(i)+(h/2),y(i,:)+(h/2)*k1)';
k3=f(t(i)+(h/2),y(i,:)+(h/2)*k2);
k4=f(t(i)+h,y(i,:)+h*k3);
y(i+1,:)=y(i,:)+h*((k1+2*k2+2*k3+k4)/6);%RK4 iteration
t(i+1,1)=t(i)+h;
end
figure(1)
plot(t,y(:,1),t,y(:,2))
legend('y1','y2')
figure(2)
plot(y(:,1),y(:,2))
xlabel('y1');
ylabel('y2');
  2 件のコメント
Riley McFerran
Riley McFerran 2022 年 4 月 25 日
My issue is on line 32, or: y(i+1,:)=y(i,:)+h*((k1+2*k2+2*k3+k4)/6);%RK4 iteration
I'm not sure how to fix it
Riley McFerran
Riley McFerran 2022 年 4 月 25 日
This is the code I typed into the command window. After entering the second line, is when I get the error
pp=@(t,x)[.23*x(1)-.0133*x(1)*x(2); -.4*x(2)+.0004*x(1)*x(2)]
[t,x]=rk4sys(pp,[0,65],[550,20],0.3)

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

回答 (1 件)

VBBV
VBBV 2022 年 4 月 25 日
pp=@(t,x)[.23*x(1)-.0133*x(1)*x(2); -.4*x(2)+.0004*x(1)*x(2)]
pp = function_handle with value:
@(t,x)[.23*x(1)-.0133*x(1)*x(2);-.4*x(2)+.0004*x(1)*x(2)]
[t,y]=rk4sys(pp,[0,65],[550,20],0.3)
n = 217
t = 1×218
0 0.3000 0.6000 0.9000 1.2000 1.5000 1.8000 2.1000 2.4000 2.7000 3.0000 3.3000 3.6000 3.9000 4.2000 4.5000 4.8000 5.1000 5.4000 5.7000 6.0000 6.3000 6.6000 6.9000 7.2000 7.5000 7.8000 8.1000 8.4000 8.7000
y = 218×2
550.0000 20.0000 545.2491 18.9428 542.7727 17.9337 542.4326 16.9756 544.1123 16.0699 547.7151 15.2174 553.1617 14.4179 560.3886 13.6709 569.3460 12.9751 579.9967 12.3293
figure(1)
plot(t,y(:,1),t,y(:,2))
legend('y1','y2')
figure(2)
plot(y(:,1),y(:,2))
xlabel('y1');
ylabel('y2');
function [ t,y ] = rk4sys( f,tspan,y0,h )
%function [ t,y ] = rk4sys( f,tspan,y0,h )
%Uses RK4 method to solve the system ode y1'=p(t,y1,y2) y1(t0)=y1_0
%y2'=q(t,y1,y2) y2(t0)=y2_0
%on the interval tspan=[t0, tend] using time-step h.
%f=[p;q] and is formatted as required by ode45
%y0=[y1_0, y2_0]
%y is the output vector approximating [y1(t), y2(t)] at times stored in t.
%calculate n the number of time-steps needed to reach tend.
n=ceil((tspan(2)-tspan(1))/h)
%Store initial conditions in t and y.
y(1,1)=y0(1);
y(1,2)=y0(2);
t(1)=tspan(1);
%RK4 Iteration Loop
for i=1:n
%If on last time-step, redefine h to make final t value = tend.
if i==n
h=tspan(2)-t(i);
end
k1=f(t(i),y(i,:)).';
k2=f(t(i)+(h/2),y(i,:)+(h/2)*k1).';
k3=f(t(i)+(h/2),y(i,:)+(h/2)*k2).'; % transpose this
k4=f(t(i)+h,y(i,:)+h*k3).'; % transpose this too
y(i+1,:)=y(i,:)+h*((k1+2*k2+2*k3+k4)/6);%RK4 iteration
t(i+1)=t(i)+h;
end
% t
% y
end
You need to transpose k3 and k4 too in function rk4sys
  2 件のコメント
Riley McFerran
Riley McFerran 2022 年 4 月 25 日
Thank you so much! Just overlooked that
VBBV
VBBV 2022 年 4 月 25 日
Welcome, please accept the answer if it resolved your problem, :)

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

カテゴリ

Help Center および File ExchangeSimulink Functions についてさらに検索

タグ

製品


リリース

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by