# 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)

2 ビュー (過去 30 日間)
Riley McFerran 2022 年 4 月 25 日
コメント済み: VBBV 2022 年 4 月 25 日
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 件のコメント表示非表示 1 件の古いコメント
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 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 件のコメント表示非表示 1 件の古いコメント
VBBV 2022 年 4 月 25 日

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

R2022a

### Community Treasure Hunt

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

Start Hunting!