四阶变步长龙格库塔法改错求助

煜宸0922

各位大侠，在下编了个变步长的四阶龙格库塔解方程组的程序，可是老出错，请各位大侠指正，程序如下：
clc
clear
AbsTol=1e-6;
h=0.01;
p=4;
a=0;
b=500;
p(21)=2^p-1;
T=a:h:b;
while T(end)<T(500);
    [T(1),X(:,1)]=rk44(a,b,h);
    [T(2),X(:,2)]=rk44(a,b,h/2);
    if abs(X(:,1)-X(:,2))/p(21)>AbsTol;
        T(1)=T(2);
        X(:,1)=X(:,2);
        h=h/2;
        [T(2),X(:,2)]=rk44(a,b,h/2);
    end
end
[T(j+1),X(:,j+1)]=rk44(a,b,h);
j=j+1;
if abs(X(1,j)-0.1)<=10e-6;
    X(3,j+1)=-0.8*X(3,j);
end
plot(X(1,(end-1000):end),X(3,(end-1000):end),'-k')
错误信息如下
??? Subscript indices must either be real positive integers or logicals.

Error in ==> changerk44 at 20
[T(j+1),X(:,j+1)]=rk44(a,b,h);请各位大侠帮忙！不胜感激

VibrationMaster

执行到[T(j+1),X(:,j+1)]=rk44(a,b,h)这句,程序不知道 j为多少

meiyongyuandeze

在你的循环开始的时候给j赋个初值！

煜宸0922

回复 3 # meiyongyuandeze 的帖子

大侠，我问的就是在这个程序基础上改进正确的变步长龙格库塔，那个rk44程序如下：
function [t,x]=rk44(a,b,h)
x(:,1)=[0 0 0 0];
t=a:h:b;
for j=1:length(t)
k1=h*feval('fun',t(j),x(:,j));
k2=h*feval('fun',t(j)+h/2,x(:,j)+k1/2);
k3=h*feval('fun',t(j)+h/2,x(:,j)+k2/2);
k4=h*feval('fun',t(j)+h,x(:,j)+k3);
x(:,j+1)=x(:,j)+(k1+2*k2+2*k3+k4)/6;
end
务求大侠帮忙，是不是while循环中不能再调用for循环，而应该先用定步长的算一个点，再在主程序中用它循环，但是这样调用我不会弄呀，所以烦求大侠赐教，不胜感激！

kezairenjian

function x=rk44(t,x,h)
k1=h*feval('fun',t,x);
k2=h*feval('fun',t+h/2,x+k1'/2);
k3=h*feval('fun',t+h/2,x+k2'/2);
k4=h*feval('fun',t+h,x+k3');
x=x+(k1'+2*k2'+2*k3'+k4')/6;
fun函数的参数为初值x，t。
是不是while循环中不能再调用for循环，而应该先用定步长的算一个点，再在主程序中用它循环？
while循环中不能再调用for循环？关键看你要实现什么。好好看看基本的东西。

daivder

我正想研究这个呢

煜宸0922

回复 5 # kezairenjian 的帖子

这位，你说的很对，就是因为while中不能再调for，我要用的是变步长的龙格库塔法，简单说，就是对下一个点同时进行以h和h/2为步长的计算，最后对两结果进行差值比较，如果差值在误差范围内就继续用h，如果超出误差范围就用h/2进行计算。大侠，如果可以，帮忙给个程序。