Why is this giving me an error message when trying to use trapz.

Hello, I'm trying to integrate for each time step the following. But it keeps giving me an error message.

B=18.6;
t0=0;
step=0.5;
tmax=90;
npas=tmax/step+1;
U=60;
for i= 1:npas
    t0=t0+step;
    s = (U.*t0)./B;
    fun3 = @(s) (0.0081*(exp(-0.0058.*s))-0.0392.*exp(-0.0833.*s));
    y3 = fun3(s);
    q3(i) = trapz(s,y3);
end
B=18.6;
t0=0;
step=0.5;
tmax=90;
npas=tmax/step+1;
U=60;
for i= 1:npas
    t0=t0+step;
    s = (U.*t0)./B;
    fun3 = @(s) (0.0081*(exp(-0.0058.*s))-0.0392.*exp(-0.0833.*s));
    y3 = fun3(s);
    q3(i) = trapz(s,y3);
end

[Geoff Hayes]

Saruultugs - is the error

Error using permute
ORDER contains an invalid permutation index.

Error in trapz (line 43)
  y = permute(y,perm);

Error in *** (line 14)
    q3(i) = trapz(s,y3);

In the future, please include the full error message.

[Geoff Hayes]

Are your inputs to trapz correct? Should y3 be a scalar?

[Saruultugs Batzorig]

Error usin
g trapz (line 47)
Dimension argument must be a positive integer scalar within indexing range.

Error in test (line 13)
    q3(i) = trapz(s,y3);

Sorry, this is the error code, I also realise I entered the code twice in my question.

[Walter Roberson]

Geoff's hint is correct:

trapz() with two parameters can be either

trapz(X, Y)

or

trapz(X, dimension)

The way it tells the two apart is that if the second parameter is a scalar, it assumes that it is a dimension.

You should not be using trapz step by step: you should be recording all of those s and y3 values and doing a trapz() on the final vector of them.

[Walter Roberson]

fun3 = @(s) (0.0081*(exp(-0.0058.*s))-0.0392.*exp(-0.0833.*s));  %does not change with i

for i= 1:npas
    t0=t0+step;
    s(i) = (U.*t0)./B;
    y3(i) = fun3(s(i));  
end
q3 = trapz(s, y3);

Though I suspect you would prefer to use cumtrapz() instead of trapz() here.

[Saruultugs Batzorig]

It worked! Thank you so much.

Answer 1

[Allen]

A more concise aproach.

B = 18.6;
step = 0.5;
tmax = 90;
npas = tmax/step+1;
U = 60;
fun3 = @(s) (0.0081*(exp(-0.0058.*s))-0.0392.*exp(-0.0833.*s));
s = U/B*step*(1:npas/2/step);
y3 = fun3(s);
q3 = trapz(s,y3);