matlab 微分方程组的向量编码及其近似解

6za6bjd0  于 2023-01-05  发布在  Matlab
关注(0)|答案(1)|浏览(189)

我试图用欧拉方法求解一个微分方程组,首先,我将系统编码成一个向量,然后将初始条件传递给函数ode_Euler。
但是,我的尝试出现了一些错误。我收到以下错误:

>> 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;

这是我的代码:

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

下面是欧拉方法:

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

当我使用另一个方法ode_Heun时,我得到了同样的错误:

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

我可以得到一些帮助来了解我的代码问题吗?

t1rydlwq

t1rydlwq1#

y(k,:)是一个行向量,而f的返回值是一个列向量。根据广播规则,行向量和列向量之和是一个矩阵,即重复的行向量和列向量的矩阵之和。
这在向量和矩阵运算的上下文中不是很合乎逻辑,但是在处理向量的(有限)序列时是有意义的,不幸的是,这种区别在类型系统中没有实现和实施。

相关问题