I am trying to solve a system of differential equations by the Euler's method. First, I have encoded my system into a vector. Then I pass the initial conditions to the function ode_Euler.
However, there's somehting wrong about my attempt. I am getting this error:
>> nm06p03a
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.
Error in ode_Euler (line 11)
y(k 1,:)= y(k,:) h*feval(f,t(k),y(k,:));
Error in nm06p03a (line 12)
tic, [tE,xE]=ode_Euler(f,tspan,x0,N); t_Euler=toc;
This is my code so far:
clear, clf
f=@(t,x)[-x(2) x(1)*x(2); x(1)-(0.5.*x(1).^2) (0.5.*x(2).^2)]; %Encoding system of differential equations into a vector
t0=0; tf=10;
tspan=[t0 tf];
N=100;
x0s=[0.2 0]; % A matrix consisting of initial points
for iter=1:size(x0s,1)
x0=x0s(iter,:);
tic, [tE,xE]=ode_Euler(f,tspan,x0,N); t_Euler=toc;
subplot(220 iter),
plot(tE,xE,'r:')
legend('ode_ Euler')
end
Here is the Euler's method:
function [t,y]=ode_Euler(f,tspan,y0,N)
if nargin<4|N<=0, N=100; end
if nargin<3, y0=0; end
h=(tspan(2)-tspan(1))/N;
t=tspan(1) [0:N]'*h;
y(1,:)=y0(:)'; %make it a row vector
for k=1:N
y(k 1,:)= y(k,:) h*feval(f,t(k),y(k,:));
end
When I use this other method ode_Heun
, I get the same error:
function [t,y]=ode_Heun(f,tspan,y0,N)
if nargin<4|N<=0, N=100; end
if nargin<3, y0=0; end
h=(tspan(2)-tspan(1))/N; % Step-size
t=tspan(1) [0:N]'*h; % Time vector
y(1,:)=y0(:)'; % make the initial value a row vector
for k=1:N
fk= feval(f,t(k),y(k,:)); y(k 1,:)= y(k,:) h*fk; % Eq.(6.2.3)
y(k 1,:)= y(k,:) h/2*(fk feval(f,t(k 1),y(k 1,:))); % Eq.(6.2.4)
end
Can I get some help to understand the problem with my code?
CodePudding user response:
y(k,:)
is a row vector, while the return value of f
is a column vector. Per broadcasting rules the sum of a row and a column vector is a matrix as the sum of the matrices of repeated row and column vectors.
This is not very logical in the context of vector and matrix operations, but can make sense in the processing of (finite) sequences of vectors. Unfortunately that distinction is not realized and enforced in the type system.