How can I succesfully implement an adaptive Rk method without getting any error like "Subscript indices must either be real positive integers or logicals" or "Error using * Inner matrix dimensions must agree"?.

2 次查看（过去 30 天）
Alexandru Lapusneanu
评论： Alexandru Lapusneanu 2017-10-23
I'm currently working on an orbit propagator and I'm trying to implement an adaptive rk method , but it doesn't work.I've started by implementing a normal Rk4 and it's working great, but now I need a more precise method due to the fact that most of the values are not very accurate.
This is the method that works.
Now the first try of an adaptive method:
h=1;
tfin=86400;
N=(tfin/h)+1;
y = zeros(6, tfin/h);
y(:,1) = y0;
epsilon=10^(-5);
N=(tfin/h)+1;
for i=1:(N/h)
tin=h*(i-1);
H=sqrt(y(1,i)^2+y(2,i)^2+y(3,i)^2);
k_1 =h* proiectia(tin(i), y(:,i), H, m, A, y(4:6, i));
k_2 =h* proiectia(tin(i)+h/4, y(:,i)+k_1/4, H, m, A, y(4:6, i));
k_3 =h* proiectia((tin(i)+(3*h)/8), (y(:,i)+((3/32)*k_1)+((9/32)*k_2)), H, m, A, y(4:6, i));
k_4 =h* proiectia((tin(i)+(12*h)/13),(y(:,i)+((1932/2197)*k_1)-((7200/2197)*k_2)+((7296/2197)*k_3)), H, m, A, y(4:6, i));
k_5 =h* proiectia(tin(i)+h,(y(:,i)+((439/216)*k_1)-(8*k_2)+((3680/513)*k_3)-(845/4104)*k_4), H, m, A, y(4:6, i));
k_6 =h* proiectia(tin(i)+h/2,(y(:,i)-((8/27)*k_1)+(2*k_2)-((3544/2565)*k_3)+(1859/4104)*k_4-((11/40)*k_5)), H, m, A, y(4:6, i));
w1 = y(:,i) + ((25/216)*k_1)+((1408/2565)*k_3)+((2197/4104)*k_4)-((1/5)*k_5);
w2 = y(:,i) + ((15/135)*k_1)+((6656/12825)*k_3)+((28561/56430)*k_4)-((9/50)*k_5)+((2/55)*k_6);
R=abs(w1-w2);
delta=0.84*(epsilon/R).^(1/4);
if R<=epsilon
tin=tin+h;
y(:,i)=w1;
h=delta*h;
else
h=delta*h;
end
end
The error was:
Error using *
Inner matrix dimensions must agree.
Can someone help me understand how to solve this problem?

请先登录，再回答此问题。

回答（2 个）

Alexandru Lapusneanu
编辑：Walter Roberson 2017-10-23
%
function yprim=proiectia (t,y,H,m,A,vit)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%Acceleratia datorata Drag-ului Atmosferic
%Parametri
ro=1.225*exp(-sqrt((y(1)^2+y(2)^2+y(3)^2)^2) /8); % Densitatea
CD=2.2; % Drag coeficient
ad=(-1/2)*ro*((CD*A)/m)*norm(vit)^2*(vit/norm(vit));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Durata de rulare in functie de Ziua Juliana
%Parametri
an=2019;
luna=3;
zi=29;
ora=9;
min=20;
sec=52.9996;
sec = sec + rem(t, 60);
min = rem(min + floor(t/60), 60);
ora = rem(ora + floor(t/3600), 24);
zi = zi + floor(t/24/3600);
%Determinare zi Juliana
Jo = (367*an) - fix((7*(an+fix((luna+9)/12)))/4) + fix(275*luna/9) + zi + 1721013.5;
JD = Jo + (((((sec/60)+min)/60)+ora)/24);
%Parametri gravitationali
miu= 398900.5*10^9; % miu Pamant
%Ecuatii folosite pentru propagator
yprim=zeros(6,1);
yprim(1,1)=y(4);
yprim(2,1)=y(5);
yprim(3,1)=y(6);
yprim(4,1)=double(-miu*y(1)/(((y(1)^2)+y(2)^2+y(3)^2)^(3/2)))+ad(1);
yprim(5,1)=double(-miu*y(2)/(((y(1)^2)+y(2)^2+y(3)^2)^(3/2)))+ad(2);
yprim(6,1)=double(-miu*y(3)/(((y(1)^2)+y(2)^2+y(3)^2)^(3/2)))+ad(3);
end
Walter Roberson
Walter Roberson 2017-10-23
Your h starts off scalar.
The line
w1 = y(:,i) + ((25/216)*k_1)+((1408/2565)*k_3)+((2197/4104)*k_4)-((1/5)*k_5);
is clearly intended to work with a vector, and so w1 will give a vector result. The same is true for the next line for w2.
R=abs(w1-w2);
is subtracting one vector from another, so you are going to get a vector result.
We then encounter
delta=0.84*(epsilon/R).^(1/4);
epsilon is a scalar, so this is a scalar with Matrix Right Divide by a vector. It so happens that the orientation of R (column vector) is acceptable for Matrix Right Divide, so this gives a vector result (a row vector). If we hyphothesize that Right Divide was wanted instead of Matrix Right Divide, 0.84*(epsilon./R).^(1/4) then the result of that hypothetical operation would be a column vector. So delta is either going to be a row vector (current code) or column vector (hypothetical code.)
You then have
if R<=epsilon
we established above that R is a vector, so remember that this code is equivalent to
if all(R<=epsilon)
that is, the body of the "if" is done only if all of the R values in the vector are less than epsilon. Perhaps that really is the desired behaviour: if it is, then please code the all() yourself to make it easier for other people to read.
Both the true and the false branches of the "if" end with
h=delta*h;
We established above that delta is either a row vector (current code) or a column vector (hypothetical code). Either way, h ends up being a vector. The length of the vector will, in either case, end up being the same as the length returned by proiectia
Next loop iteration, we encounter
k_1 =h* proiectia(tin(i), y(:,i), H, m, A, y(4:6, i));
but now h is a vector. The routine proiectia returns a column vector. Can the * operation succeed?
If h is a row vector (current code) then that is 1 x N (from h) * N x 1 (from proiectia) which is going to give a scalar. You get scalars for the k_* but when you get down to the w1 and w2, the y(:,i) force a vector return, so again R is going to be a column vector. Working with the hypothesis that we continue using the current code epsilon/R rather than the hypothetical code epsilon./R then the delta that is produced will again be a row vector. Then in either branch of the "if" you delta*h . delta is a row vector, h is (under the current code) a row vector, so you have 1 x N * 1 x N and that gives you a Matrix Times error because the matrix dimensions do not agree.
If, though, h is a column vector (under the hypothesis that you want epsilon./R) then the k_1 line is N x 1 (hypothetical h) * N x 1 (from proiectia), and that gives you a Matrix Times error because the matrix dimensions do not agree.
Therefore, regardless of whether the code should be epsilon/R or elpsilon./R, the rest of the code is broken as to how it uses the * operator.
You should be stepping through the code operator by operator asking yourself what size you expect a variable to be, and whether you expect it might change size later, and what orientation you want it to be, and whether you intend to do vector "if" tests, and whether you are knowingly using Matrix Right Divide or whether you want (element-by-element) Right Divide instead.
When writing code involving vectors, I find it safer to code the dot forms of operations except where I deliberately intend to use a matrix algebraic operation. It doesn't hurt for scalar .* or ./ or .^ and doing this consistently catches mistakes in a consistent way. Especially for the / operator, which sometimes accidentally matches on size, silently giving results that mean something entirely different than what I was thinking of.

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by