How do I calculate an integral depending on solutions of an ODE system solved by Euler's method?

Hi, I want to calculate an integral which depends on the solutions of an ODE system. This ODE system is composed by 8 equations and I solved it using Euler's method. This solutions were saved in a vector of length "vectorsize". Here's the code:

ttotal=30; deltat=0.1; vectorsize=300;

deltas=0.1; vectorsizes = 300;

%VECTORS t = zeros(1,vectorsize); Sm = zeros(1,vectorsize); Sf = zeros(1,vectorsize); If = zeros(1,vectorsize); Im = zeros(1,vectorsize); Sd =zeros(1,vectorsize); Id = zeros(1,vectorsize); Sc = zeros(1,vectorsize); Ic =zeros(1,vectorsize); Tm = zeros(1,vectorsizes); Tf = zeros(1,vectorsizes);

%INITIAL CONDITIONS Sm(1)=33631; Im(1)=5248; Sf(1)=31674; If(1)=10446; Sd(1)=445; Id(1)=0; Sc(1)=444; Ic(1)=0;

Am=173; Af=187; miu=6.67*10^(-4); miu1=2.38*10^(-3); v=6.1*10^(-3);

Bf=0.0785; Bd0=0.0047; Bc0=0.0085; Bd=0.0219; Bc=0.0396; Bm=0.0341;

%EULER'S METHOD for i=2:vectorsize

    dSm = Am-((Bm*If(i-1)*Sm(i-1))/(Sf(i-1)+If(i-1)))-((miu+miu1)*Sm(i-1));
    dIm = ((Bm*If(i-1)*Sm(i-1))/(Sf(i-1)+If(i-1)))-((miu+miu1+v)*Im(i-1));
    dSf = Af-((Bf*Im(i-1)*Sf(i-1))/(Sm(i-1)+Im(i-1)))-((miu+miu1)*Sf(i-1));
    dIf = ((Bf*Im(i-1)*Sf(i-1))/(Sm(i-1)+Im(i-1)))-((miu+miu1+v)*If(i-1));

    Sm(i) = Sm(i-1)+dSm*deltat;
    Im(i) = Im(i-1)+dIm*deltat;
    Sf(i)= Sf(i-1)+dSf*deltat;
    If(i) = If(i-1)+dIf*deltat; 

    dSd = -((Bd*Im(i-1)*Sd(i-1))/(Sm(i-1)+Im(i-1)));
    dId = ((Bd*Im(i-1)*Sd(i-1))/(Sm(i-1)+Im(i-1)));
    dSc = -((Bc*Im(i-1)*Sc(i-1))/(Sm(i-1)+Im(i-1)));
    dIc = ((Bc*Im(i-1)*Sc(i-1))/(Sm(i-1)+Im(i-1)));

    Sd(i) = Sd(i-1)+dSd*deltat;
    Id(i) = Id(i-1)+dId*deltat;
    Sc(i) = Sc(i-1)+dSc*deltat;
    Ic(i) = Ic(i-1)+dIc*deltat;

    t(i)=t(i-1)+deltat;

end

And then, I wanted to calculate the integral of functions Tm and Tf (baring in mind that I solved If(s),Sm(s), Sf(s) and Im(s) with Euler's method and saved the solutions in their respective vectors):

Tm= @(s) (Bm*If(s)*Sm(s))/(Sf(s)+If(s))*deltas; Tf =@(s) (Bf*Im(s)*Sf(s))/(Sm(s)+Im(s))*deltas;

integralTm=integral(Tm,0,ttotal); integralTf=integral(Tf,0,ttotal);

------------------------------------------------------------

And then this ERROR appears:

Subscript indices must either be real positive integers or logicals.

Error in @(s)(Bm*If(s)*Sm(s))/(Sf(s)+If(s))

Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);

Error in integralCalc/vadapt (line 133) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 76) [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 89) Q = integralCalc(fun,a,b,opstruct);

Error in msm2 (line 182) integralTm=integral(Tm,0,ttotal);

--------------------------

Thank you in advance, hope you can help me out! P.S.: If necessary, I can exchange Euler's method to ode45, but only if is really necessary for this integrals to be solved.

Answer 1

[Star Strider]

The essence: ‘And then, I wanted to calculate the integral of functions Tm and Tf (bearing in mind that I solved If(s),Sm(s), Sf(s) and Im(s) with Euler's method and saved the solutions in their respective vectors):’

That also means that Tm and Tf are vectors rather than functions. Define them as:

Tm =  (Bm.*If.*Sm./(Sf+If).*deltas; 
Tf =  (Bf.*Im.*Sf./(Sm+Im).*deltas;

NOTE: For this to work, the vectors all have to be the same size.

Then use trapz to integrate them:

integralTm = trapz(Tm); 
integralTf = trapz(Tf);

Does that do what you want?

(BTW, I suggest the MATLAB ODE solvers. They’re much more accurate than Euler’s method.)