Using of cumprod fuction

I want to calculate the following formula
So I wrote out the following code.
syms Nb Nv K n Tb Hb t Sb F Sv Hv...
D P Cj j r q
Tb=1.2;
Hb=[2.8 2.2 1.7 1.2 0.8 0.6 0 0];
Sb=[25 28 32 36 35 38 0 0];
F=[70 50 48 55 80 65 0 0];
Sv=[0 280 200 320 360 380 420 0];
Hv=[0 2.5 1.9 1.4 1 0.7 0.5 0];
D=[180000 200000 240000 300000 350000 400000 0 0];
P=[0 200000 240000 300000 350000 400000 450000 0];
q=[1 1 1 1 1 1 0 0];
u=[0.009 0.0064 0.00054 0.0072 0.01 0.008 0 0];
t=[0.015 0.01 0.009 0.012 0.016 0.014 0 0];
w=[0.005 0.004 0.003 0.004 0.006 0.004 0 0];
Se=[29 32 36 35 38 0 0];
n=[1 0.9 0.833333333333333 0.8 0.857142857142857 0.875 0 0];
p=[3.2 2.6 2.1 1.6 1.2 0.9 0];
%
Nb=[1 1 1 1 1 1 0 0];
Nv=[0 8 1 1 1 1 1 0];
K=[0 10 11 9 6 4 4 0];
m=[0.0096 0.0079 0.0072 0.0092 0.0118 0.0105 0 0];
%
a = numel(Nv);
Z = sum(n(a)*Tb.*Hb(a).*D(a).*(m(a)-u(a))+ ...
(1./(cumprod(Nv(2:a).*K(a+1))).*(n(a).*Tb)^2* ...
((D(a)/2).*Hv(a+1)-(2*Hv(a+1)*(1-(D(a)/(P(a+1)))))))+ ...
(cumprod(Nv(2:a)).*K(a+1).*((Se(a).*((t(a)-m(a))^3.*(t(a)+3*m(a)-4*w(a)) ...
./(t(a)-w(a))^4))+((Hb(a)+p(a).*D(a)).*((t(a)-m(a))^5.*(2*t(a)+3*m(a)-5*w(a)) ...
./5*(t(a)-w(a))^4))+F(a)))+cumprod(Nv(2:a)).*Sv(a+1)+...
(1./cumprod(Nv(2:a))*(n(a)*Tb)^2).*((Hv(a+1)*D(a)/2).*(1-(D(a).*(p(a+1)))))...
+cumprod(Nv(1:a)).*Nb(a).*Sb(a));
format long G
Z
It display error using below.
Index exceeds the number of array elements. Index must not exceed 8.
Error in paper (line 35)
(1./(cumprod(Nv(2:a).*K(a+1))).*(n(a).*Tb)^2* ...
I think my use of cumprod is wrong. There is only 7 in the sequence, but after defining a, I have to change it to 8.
And the answer to z will be 81175.
Please help me fix the code so that I can calculate the correct result, thank you for taking the time to answer.

 Accepted Answer

Torsten
Torsten on 2 Apr 2024
Edited: Torsten on 2 Apr 2024
Tb=1.2;
Hb=[2.8 2.2 1.7 1.2 0.8 0.6 0 0];
Sb=[25 28 32 36 35 38 0 0];
F=[70 50 48 55 80 65 0 0];
Sv=[0 280 200 320 360 380 420 0];
Hv=[0 2.5 1.9 1.4 1 0.7 0.5 0];
D=[180000 200000 240000 300000 350000 400000 0 0];
P=[0 200000 240000 300000 350000 400000 450000 0];
q=[1 1 1 1 1 1 0 0];
u=[0.009 0.0064 0.00054 0.0072 0.01 0.008 0 0];
t=[0.015 0.01 0.009 0.012 0.016 0.014 0 0];
w=[0.005 0.004 0.003 0.004 0.006 0.004 0 0];
Se=[29 32 36 35 38 0 0];
n=[1 0.9 0.833333333333333 0.8 0.857142857142857 0.875 0 0];
p=[3.2 2.6 2.1 1.6 1.2 0.9 0];
%
Nb=[1 1 1 1 1 1 0 0];
Nv=[0 8 1 1 1 1 1 0];
K=[0 10 11 9 6 4 4 0];
m=[0.0096 0.0079 0.0072 0.0092 0.0118 0.0105 0 0];
%
s = 0;
for i = 1:6
s = s + n(i)*Tb*Hb(i)*D(i)*(m(i)-u(i))+ ...
(n(i)*Tb)^2/(prod(Nv(2:i+1))*K(i+1))* ...
(D(i)/2*(Hv(i+1)-2*Hv(i+1)*(1-D(i)/P(i+1))))+ ...
prod(Nv(2:i+1))*K(i+1)*(Se(i)*(t(i)-m(i))^3*(t(i)+3*m(i)-4*w(i)) ...
/(t(i)-w(i))^4 +(Hb(i)+p(i))*D(i)*(t(i)-m(i))^5*(2*t(i)+3*m(i)-5*w(i)) ...
/(5*(t(i)-w(i))^4)+F(i))+prod(Nv(2:i+1))*Sv(i+1)+...
(n(i)*Tb)^2/prod(Nv(2:i+1))*(Hv(i+1)*D(i)/2*(1-D(i)/P(i+1)))...
+ prod(Nv(2:i+1))*Nb(i)*Sb(i)
end
s = 1.6831e+04
s = 3.0432e+04
s = 4.5243e+04
s = 5.6491e+04
s = 6.8221e+04
s = 7.8814e+04
idx = 1:6;
Nvc = cumprod(Nv(idx+1));
s = sum(n(idx).*Tb.*Hb(idx).*D(idx).*(m(idx)-u(idx))+...
(n(idx)*Tb).^2./(Nvc.*K(idx+1)).* ...
(D(idx)/2.*(Hv(idx+1)-2*Hv(idx+1).*(1-D(idx)./P(idx+1))))+ ...
Nvc.*K(idx+1).*(Se(idx).*(t(idx)-m(idx)).^3.*(t(idx)+3*m(idx)-4*w(idx)) ...
./(t(idx)-w(idx)).^4 +(Hb(idx)+p(idx)).*D(idx).*(t(idx)-m(idx)).^5.*(2*t(idx)+3*m(idx)-5*w(idx)) ...
./(5*(t(idx)-w(idx)).^4)+F(idx))+Nvc.*Sv(idx+1)+...
(n(idx)*Tb).^2./Nvc.*(Hv(idx+1).*D(idx)/2.*(1-D(idx)./P(idx+1)))...
+ Nvc.*Nb(idx).*Sb(idx))
s = 7.8814e+04

1 Comment

Although it is different from the answer z, thank you for your answer. I will look again to see if there are any typos.

Sign in to comment.

More Answers (0)

Products

Release

R2023a

Tags

Asked:

on 2 Apr 2024

Commented:

on 3 Apr 2024

Community Treasure Hunt

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

Start Hunting!