我正在尝试实现中点公式y[n+1/2] = y[n] + h/2 f (x[n], y[n])
和y[n+1] = y[n] + h *f (x[n] + h/2, y[n + 1/2])
因此它使用中点方法求解ODE。
我的职能是
function [ x, y ] = Midpoint_ODE ( f, xRange, yInitial, numSteps )
% f = name of file with function
% xrange Interval
% x(1) first meaning of x
% x(2) second meaning of x
x=zeros(numSteps+1,1);
x(1) = xRange(1);
h = ( xRange(2) - xRange(1) ) / numSteps; % calculated step size
y(1,:) = transpose(yInitial);
for n = 1 : numSteps
y(n+0.5,:)= (y(n) + (h * 0.5)*(transpose(feval( f, x(n), y(n)))));
y(n+1,:) = y(n,:) + h * transpose(feval(f, x(n)+ (h/2), y(n+0.5,:))); %evaluating the function
end
但我收到一个错误:
**Index in position 1 is invalid. Array indices must be positive integers or logical values.
Error in Midpoint_ODE (line 11)Index in position 1 is invalid. Array indices must be positive integers or logical values.
Error in Midpoint_ODE (line 11)**
我检查了几次,但如果我错过了一些合乎逻辑的内容,就无法理解问题所在。
您不需要保留半步值。因此,最简单的是不在值列表中]
for n = 1 : numSteps yhalfstep = (y(n,:) + (h * 0.5)*(transpose(feval( f, x(n), y(n,:))))); y(n+1,:) = y(n,:) + h * transpose(feval( x(n)+ (h/2), yhalfstep)); end
还请记住,在matlab和类似的软件中,对多维数组的单索引访问会返回展平数组的元素(列优先)。也就是说,在
a=[ 1,2;3,4;5,6]
中,您从a(3)
中获得数字5
作为第一列中的第三元素,而a(3,:)
给出第三行[5,6]
。