MATLAB: How to calculate an integral depending on solutions of an ODE system solved by Euler’s method

integrationnumerical integrationode45vectors

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.

Best Answer

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.)
Related Question