Octave ODE solver: inconsistent size for state and derivative vectors

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Octave ODE solver: inconsistent size for state and derivative vectors

Fabcap77
I have the code to solve a problem of non-stationary heat transfer (no, this
is not homework, but the code is taken from a textbook), which requires
solving a set of ODEs with the Method of Lines, and goes as follows (I added
some diagnostic instructions):

%Problem P6_08A
clear, clc, format short g, format compact
 c_time = clock();
 r_time = r_time =
strjoin({num2str(c_time(4)),num2str(c_time(5)),num2str(c_time(6))}, ':');
 disp('Run starts at'); disp(r_time);
tspan = [0 10]; % Range for the independent variable
y0 = [100.; 100.; 100.; 100.; 100.; 100.; 100.; 100.; 100.]; % Initial
values for the dependent variables
%y0 = 100*ones(1000,1);
%- - - - - - - - - - - - - - - - - - - - - -
function dYfuncvecdt = ODEfun(Yfuncvec, t);
 disp('Inside ODEfun');
%Yfuncvec = 100*ones(9,1);
T2 = Yfuncvec(1);
T3 = Yfuncvec(2);
T4 = Yfuncvec(3);
T5 = Yfuncvec(4);
T6 = Yfuncvec(5);
T7 = Yfuncvec(6);
T8 = Yfuncvec(7);
T9 = Yfuncvec(8);
T10 = Yfuncvec(9);
 more off;
 disp('Yfuncvec = '); disp(Yfuncvec);
alpha = .00002;
deltax = .1;
T1 = 0;
T11 = (4 * T10 - T9) / 3;
dT2dt = alpha / (deltax ^ 2) * (T3 - (2 * T2) + T1);
dT3dt = alpha / (deltax ^ 2) * (T4 - (2 * T3) + T2);
dT4dt = alpha / (deltax ^ 2) * (T5 - (2 * T4) + T3);
dT5dt = alpha / (deltax ^ 2) * (T6 - (2 * T5) + T4);
dT6dt = alpha / (deltax ^ 2) * (T7 - (2 * T6) + T5);
dT7dt = alpha / (deltax ^ 2) * (T8 - (2 * T7) + T6);
dT8dt = alpha / (deltax ^ 2) * (T9 - (2 * T8) + T7);
dT9dt = alpha / (deltax ^ 2) * (T10 - (2 * T9) + T8);
dT10dt = alpha / (deltax ^ 2) * (T11 - (2 * T10) + T9);
dYfuncvecdt = [dT2dt; dT3dt; dT4dt; dT5dt; dT6dt; dT7dt; dT8dt; dT9dt;
dT10dt];
end
%Solution of ODE set
 y = lsode("ODEfun", y0,tspan);
%
disp([y0 ODEfun(tspan(1),y0)]);
disp(' Variable values at the initial point ');
disp([' t    = ' num2str(tspan(1))]);
disp('           y                  dy/dt         ');
for i=1:size(y,2)
    disp([' Solution for dependent variable y' int2str(i)]);
    disp(['              t                  y' int2str(i)]);
    disp([t y(:,i)]);
    plot(t,y(:,i));
    title([' Plot of dependent variable y' int2str(i)]);
    xlabel(' Independent variable (t)');
    ylabel([' Dependent variable y' int2str(i)]);
    pause
end
%EOF

Trying to run the code returns an error:

error: ODEfun: A(I): index out of bounds; value 2 out of bound 1
error: called from:
error:   ODEfun at line 14, column 4
error: evaluating argument list element number 1
error:  $path/lines01_diff_eq.m at line 42, column 1

I posted this also on other forums, but nobody could help. I hope this is
the right place for asking.



--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: Octave ODE solver: inconsistent size for state and derivative vectors

Colin Macdonald-2
On 2018-04-13 08:48 PM, Fabcap77 wrote:
> disp([y0 ODEfun(tspan(1),y0)]);

Should this be "ODEfun(y0, tspan(1))"?



Reply | Threaded
Open this post in threaded view
|

Re: Octave ODE solver: inconsistent size for state and derivative vectors

Fabcap77
Oh yes, that line escaped my attention. I corrected it, and after some more
tweaking the code can be executed without errors. Whether the results make
sense phisically, it will take some more investigations to tell.

Thank for the help.



--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: Octave ODE solver: inconsistent size for state and derivative vectors

Przemek Klosowski-7
In reply to this post by Fabcap77
On 04/13/2018 11:48 PM, Fabcap77 wrote:
> error: ODEfun: A(I): index out of bounds; value 2 out of bound 1
> error: called from:
> error:   ODEfun at line 14, column 4
> error: evaluating argument list element number 1
> error:  $path/lines01_diff_eq.m at line 42, column 1

Well, you call disp([y0 ODEfun(tspan(1),y0)]);

and tspan(1) is a number, while ODEfun wants an array. I didn't try to
figure out what is the intent of that statement, but you have to give it
a 9-element array as the first argument.