Symbolic integration loop: how to make less computationally intense?

I have a piece of code with three end equations relying heavily on symbolic integration of previous equations.

I want to be able to graph the three equations over an interval, let's call it theta or T from 0 to 1, for different parameter permutations. Specifically, I want to do this for different values of a parameter, let's call it n. For n=2, the code runs fine, but for higher levels of n, Matlab gets stuck as "busy" trying to compute one of the equations.

I decided to try to loop it for specific values of T (values shown as vector in code below), and it worked fine on my old computer, but that computer crashed and now the code is getting stuck as busy again on the same equation with a different computer. The loop I am using is:

                        %% loop_over_myScript
                        vector = [0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1];
                        out_vector = [];
                        out_vector2 = [];
                        out_vector3 = [];
                        for vec = vector
                            T = vec;
                            TestScript
                            % ignore the warning about growing in the loop
                            out_vector = [ out_vector, Delta_w ];
                            out_vector2 = [ out_vector2, Delta_npi ]; 
                            out_vector3 = [ out_vector3, Delta_u ];
                        end

Being looped, the "TestScript" which has the three equations at the end (mentioned in loop above as Delta_w, Delta_npi and Delta_u) is:

                %clear all
                %syms s e v c T n p p_x pmx Epx Eminpx
                s=0;
                e=0;
                c=0.2;
                n=20;
                %assume(p_x>=0 & pmx>=p_x);
                %assume(Epx>=0 & 0.5>=Epx & Eminpx>=0 & 0.5>=Eminpx);
                Vs=int(v,v,s-e,1)+int(s-e,v,0,s-e);
                pmo=Vs;
                npio=(1-T)*(Vs-c);
                p_o=(npio)/(1-T+(T*n))+c;
                wo=Vs-c;
                pmx=0.5*(1+s+e);
                npix=(1-T)*((pmx-s)*(1-(pmx-e))-(c-s));
                p_x=(e*((T*n) - T + 1)^(1/2) + s*((T*n) - T + 1)^(1/2) + ((T*n) - T + 1)^(1/2) - T^(1/2)*n^(1/2)*(e^2 - 2*e*s + 2*e + s^2 + 2*s - 4*c + 1)^(1/2))/(2*(T*n - T + 1)^(1/2));
                Gxp=1-(((1-T)*((pmx-s)*(1-(pmx-e))-(p-s)*(1-(p-e))))/(T*n*((p-s)*(1-(p-e))-(c-s))))^(1/(n-1));
                Epx=p_x+int(1-Gxp,p,p_x,pmx);
                Eminpx=p_x+int((1-Gxp)^n,p,p_x,pmx);
                wNSx=int(v,v,Epx-e,1)+int(s-e,v,0,Epx-e)-c;
                wSx=int(v,v,Eminpx-e,1)+int(s-e,v,0,Eminpx-e)-c;
                wx=((1-T)*wNSx)+(T*wSx);
                Delta_npi=npix-npio;
                Delta_w=vpa(wx-wo)
                Delta_u=Delta_w-Delta_npi;

The code is busy for hours and stuck calculating Epx and wNSx, and I think it is an integration problem, but not sure how to get around it. I did also try CoCalc online, and it said the "Romberg method" of integration couldn't find a numerical integral for Epx, if that helps.

Does anyone have any suggestions for alternative coding or approach that would either be less computationally intense (so my computer could run this as a loop) or that would resolve the integration problem for T from 0 to 1?

Thanks!

[Torsten]

To make code faster, don't use symbolic computations: use "integral" instead of "int".

[D Bird]

Took some rewriting, but got it to work. Very fast now. Thank you very much!

Below is the revised code as well as one more small question I hope you or someone can also help with that I ran into with the revised code.

s=0;
e=0;
c=0.2;
n=100;
T=0.5
%T=[0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1]
%assume(p_x>=0 & pmx>=p_x);
%assume(Epx>=0 & 0.5>=Epx & Eminpx>=0 & 0.5>=Eminpx);
fun=@(v) v;
Vs=integral(fun,0,1);
%For s-e integral
% Vs_int1=integral(fun,s-e,1); 
% fun1=@(v) s-e
% Vs_int2=integral(fun1,0,s-e)
% Vs=Vs_int1+Vs_int2
pmo=Vs;
npio=(1-T)*(Vs-c);
p_o=(npio)/(1-T+(T*n))+c;
%Gop=1-(((1-T)*(Vs-p))/(T*n*(p-c)))^(1/(n-1))
%Epo=p_o+int(1-Gop,p,p_o,pmo)
%Eminpo=p_o+int((1-Gop)^n,p,p_o,pmo)
wo=Vs-c;
pmx=0.5*(1+s+e);
npix=(1-T)*((pmx-s)*(1-(pmx-e))-(c-s));
p_x=(e*((T*n) - T + 1)^(1/2) + s*((T*n) - T + 1)^(1/2) + ((T*n) - T + 1)^(1/2) - T^(1/2)*n^(1/2)*(e^2 - 2*e*s + 2*e + s^2 + 2*s - 4*c + 1)^(1/2))/(2*(T*n - T + 1)^(1/2));
Gxp=@(p) (1-(1-(((1-T).*((pmx-s).*(1-(pmx-e))-(p-s).*(1-(p-e))))./(T.*n.*((p-s).*(1-(p-e))-(c-s)))).^(1./(n-1)))) %1-Gxp
%fun2=@(p) 1-Gxp
Epx_int=integral(Gxp,p_x,pmx);
Epx=p_x+Epx_int;
fun3=@(p) ((1-(1-(((1-T).*((pmx-s).*(1-(pmx-e))-(p-s).*(1-(p-e))))./(T.*n.*((p-s).*(1-(p-e))-(c-s)))).^(1./(n-1))))).^n %(1-Gxp)^n
Eminpx_int=integral(fun3,p_x,pmx);
Eminpx=p_x+Eminpx_int;
wNSx=integral(fun,Epx-e,1)-c;
%For s-e integral
% wNSx_1=integral(fun,Epx-e,1)
% wNSx_2=integral(fun1,0,Epx-e)
% wNSx=wNSx_1+wNSx_2-c;
wSx=integral(fun,Eminpx-e,1)-c;
%For s-e integral
% wSx_1=integral(fun,Eminpx-e,1)
% wSx_2=integral(fun1,0,Eminpx-e)
% wSx=wSx_1+wSx_2-c;
wx=((1-T)*wNSx)+(T*wSx);
Delta_npi=npix-npio;
Delta_w=vpa(wx-wo)
Delta_u=Delta_w-Delta_npi;

The minor question: note the commented out bits above about one of my integrals (the %For s-e integral lines that under current parameter assumptions equals zero (s-e) so I was able to just delete for now). However, if I want to include this integral, when I try to run it with this function (fun1), I get the error:

Error using integralCalc/finalInputChecks
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to
true.

Error in integralCalc/iterateScalarValued (line 315)
finalInputChecks(x,fx);

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

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

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

Error in TestScript (line 16)
Vs_int2=integral(fun1,0,s-e)

I had this error before and managed to fix it, but I can't figure out where I went wrong in the coding on this one. Any suggestions?

[Star Strider]

Setting 'ArrayValued',1 where indicated —

s=0;
e=0;
c=0.2;
n=100;
T=0.5

T = 0.5000

%T=[0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1]
%assume(p_x>=0 & pmx>=p_x);
%assume(Epx>=0 & 0.5>=Epx & Eminpx>=0 & 0.5>=Eminpx);
fun=@(v) v;
Vs=integral(fun,0,1, 'ArrayValued',1)

Vs = 0.5000

%For s-e integral
% Vs_int1=integral(fun,s-e,1); 
% fun1=@(v) s-e
% Vs_int2=integral(fun1,0,s-e)
% Vs=Vs_int1+Vs_int2
pmo=Vs;
npio=(1-T)*(Vs-c);
p_o=(npio)/(1-T+(T*n))+c;
%Gop=1-(((1-T)*(Vs-p))/(T*n*(p-c)))^(1/(n-1))
%Epo=p_o+int(1-Gop,p,p_o,pmo)
%Eminpo=p_o+int((1-Gop)^n,p,p_o,pmo)
wo=Vs-c;
pmx=0.5*(1+s+e);
npix=(1-T)*((pmx-s)*(1-(pmx-e))-(c-s));
p_x=(e*((T*n) - T + 1)^(1/2) + s*((T*n) - T + 1)^(1/2) + ((T*n) - T + 1)^(1/2) - T^(1/2)*n^(1/2)*(e^2 - 2*e*s + 2*e + s^2 + 2*s - 4*c + 1)^(1/2))/(2*(T*n - T + 1)^(1/2));
Gxp=@(p) (1-(1-(((1-T).*((pmx-s).*(1-(pmx-e))-(p-s).*(1-(p-e))))./(T.*n.*((p-s).*(1-(p-e))-(c-s)))).^(1./(n-1)))) %1-Gxp

Gxp = function_handle with value:
    @(p)(1-(1-(((1-T).*((pmx-s).*(1-(pmx-e))-(p-s).*(1-(p-e))))./(T.*n.*((p-s).*(1-(p-e))-(c-s)))).^(1./(n-1))))

%fun2=@(p) 1-Gxp
Epx_int=integral(Gxp,p_x,pmx, 'ArrayValued',1)

Epx_int = 0.2094

Epx=p_x+Epx_int;
fun3=@(p) ((1-(1-(((1-T).*((pmx-s).*(1-(pmx-e))-(p-s).*(1-(p-e))))./(T.*n.*((p-s).*(1-(p-e))-(c-s)))).^(1./(n-1))))).^n %(1-Gxp)^n

fun3 = function_handle with value:
    @(p)((1-(1-(((1-T).*((pmx-s).*(1-(pmx-e))-(p-s).*(1-(p-e))))./(T.*n.*((p-s).*(1-(p-e))-(c-s)))).^(1./(n-1))))).^n

Eminpx_int=integral(fun3,p_x,pmx, 'ArrayValued',1)

Eminpx_int = 0.0044

Eminpx=p_x+Eminpx_int;
wNSx=integral(fun,Epx-e,1, 'ArrayValued',1)-c

wNSx = 0.1815

%For s-e integral
% wNSx_1=integral(fun,Epx-e,1)
% wNSx_2=integral(fun1,0,Epx-e)
% wNSx=wNSx_1+wNSx_2-c;
wSx=integral(fun,Eminpx-e,1, 'ArrayValued',1)-c

wSx = 0.2603

%For s-e integral
% wSx_1=integral(fun,Eminpx-e,1)
% wSx_2=integral(fun1,0,Eminpx-e)
% wSx=wSx_1+wSx_2-c;
wx=((1-T)*wNSx)+(T*wSx);
Delta_npi=npix-npio;
Delta_w=vpa(wx-wo)

Delta_w = 

Delta_u=Delta_w-Delta_npi;

.

[D Bird]

Perfect! Thank you. That worked. As you could maybe tell, I am a novice coder, so thanks for the help :)

[Star Strider]

Our pleasure!