How to find the upper limit of a complex integral

I'm trying to solve this complex integral using matlab where Fa0,Fb0,Fc0,Fd0,v0 are all constants but k1 and k2 change in a for loop by pulling from a variable set. Where I want to find x and it lies between 0 and 1. The graph I want to get is this. Some clairifcation Tvals are the temperature which k1 and k2 are calculated and the x-axis of my plot. Second k1 and k2 in the list of 106 are used together. So at the first entry of Tval is 450 the first entry of k1 and k2 would be used in equation together.

load("Tvals.mat")
kvals=load("kvals.mat")
P=506625
R=8.314
Tin=750
v0=0.019635
yca0=0.1
ycb0=0.25
ycc0=0.10
ycd0=0.05
ca0=(P/(R*Tin))*yca0
cb0=(P/(R*Tin))*ycb0
cc0=(P/(R*Tin))*ycc0
cd0=(P/(R*Tin))*ycd0
Fa0=ca0*v0
Fb0=cb0*v0
Fc0=cc0*v0
Fd0=cd0*v0
xvals=[]
for i=1:length(Tvals)
k1=kvals.K1vals(i)
k2=kvals.K2vals(i)
func=@(x)(Fa0./(k1.*Fa0.*(1 - x).*(-Fa0.*x + Fb0)./v0.^2 - k2.*(Fa0.*x + Fc0).*(Fa0.*x + Fd0)./v0.^2))
xval=fzero(func,0);
xvals=[xvals xval];
end
plot(Tvals,xvals)
end

I've also tried to integrate it first

func=@(x)(Fa0./(k1.*Fa0.*(1 - x).*(-Fa0.*x + Fb0)./v0.^2 - k2.*(Fa0.*x + Fc0).*(Fa0.*x + Fd0)./v0.^2))
intfun=@(x)integral(func,0,x)-0.5
xval=fzero(intfun,0);

Answer 1

[Star Strider]

I’ve been working on this for a while.

This works, however at some point, fsolve finds undefined values at the intial point, and stops:

D1 = load('kvals.mat');
K1vals = D1.K1vals;
K2vals = D1.K2vals;
D2 = load('Tvals.mat');
Tvals = D2.Tvals;

P=506625;
R=8.314;
Tin=750;
v0=0.019635;
yca0=0.1;
ycb0=0.25;
ycc0=0.10;
ycd0=0.05;
ca0=(P/(R*Tin))*yca0;
cb0=(P/(R*Tin))*ycb0;
cc0=(P/(R*Tin))*ycc0;
cd0=(P/(R*Tin))*ycd0;
Fa0=ca0*v0;
Fb0=cb0*v0;
Fc0=cc0*v0;
Fd0=cd0*v0;
xval=zeros(1,numel(Tvals));
func=@(x,k1,k2) (Fa0./(k1.*Fa0.*(1 - x).*(-Fa0.*x + Fb0)./v0.^2 - k2.*(Fa0.*x + Fc0).*(Fa0.*x + Fd0)./v0.^2));
for k1 = 1:numel(K1vals)
    for k2 = 1:numel(K2vals)
        [xval(k1,k2),fval(k1,k2)] = fsolve(@(x)integral(@(x)func(x,K1vals(k1),K2vals(k2))-0.5,0,x),1000);
    end
end
[K2,K1] = ndgrid(K2vals, K1vals);
T1 = table(K1(:), K2(:), xval(:), fval(:),'VariableNames',{'K1','K2','Xval','Fval'});
save('How to find the upper limit of a complex integrall_T1.mat','T1');

This is my latest attempt, I have no idea if it’s going to work, since it takes forever for the two nested loops.

Note that it is necessary to iterate through all combinations of ‘K1’ and ‘K2’, since the integral changes depending on those values. If the loops complete, ‘T1’ gets saved as a .mat file so it is only necessary to do this once.

I have no idea what ‘Tvals’ does, since it does not appear other than in the plot. (Looping through it simply produces 106 repititions of the ‘K1’ and ‘K2’ loops.)

I’ll post this now, and if the loop completes, I’ll edit this and attach the .mat file.

[Confidential Confidential]

• Thanks for all your hard work yeah this one is a head scratcher... I find the fault tolerance on fsolver isn't high enough with the numbers being so tiny

This integral is based on a temperature set for a reactor and for each temperature a k1 and k2 value are generated which are rate constants for the reaction rate. X is the conversion of this reactor which I'm trying to find for each temperature...the integral defines all of this. It's turning out to be impossible to solve in matlab lol.

[Star Strider]

My pleasure!

It’s still iterating (after about 30 minutes — I should have timed it). I’ll post the results if it appears to complete successfully. Turns out that an initial estimate of 1000 is too high. Trying it again with 100, and timing it.

EDIT — (19 Nov 2020 at 2:30)

I didn’t realise that the ‘Tvals’, ‘K1vals’ and ‘K2vals’ needed to go together. This takes a bit over 20 seconds on my desktop:

xval=zeros(1,numel(Tvals));
func=@(x,k1,k2) (Fa0./(k1.*Fa0.*(1 - x).*(-Fa0.*x + Fb0)./v0.^2 - k2.*(Fa0.*x + Fc0).*(Fa0.*x + Fd0)./v0.^2));
tic
for k = 1:numel(Tvals)
    k1 = K1vals(k);
    k2 = K2vals(k);
    [xval(k),fval(k)] = fsolve(@(x)integral(@(x)func(x,k1,k2),0,x)-0.5,1E-5);
end
toc
T1 = table(Tvals(:), K1vals(:), K2vals(:), xval(:), fval(:),'VariableNames',{'T','K1','K2','Xval','Fval'});xval=zeros(1,numel(Tvals));
func=@(x,k1,k2) (Fa0./(k1.*Fa0.*(1 - x).*(-Fa0.*x + Fb0)./v0.^2 - k2.*(Fa0.*x + Fc0).*(Fa0.*x + Fd0)./v0.^2));
tic
for k = 1:numel(Tvals)
    k1 = K1vals(k);
    k2 = K2vals(k);
    [xval(k),fval(k)] = fsolve(@(x)integral(@(x)func(x,k1,k2),0,x)-0.5,1E-5);
end
toc
T1 = table(Tvals(:), K1vals(:), K2vals(:), xval(:), fval(:),'VariableNames',{'T','K1','K2','Xval','Fval'});

and produces ‘T1’ (here sampled every 5 rows):

     T          K1            K2           Xval          Fval   
    ____    __________    __________    __________    __________

     450    8.7311e-12    2.1726e-14    4.5138e-09    1.1102e-16
     500    9.3862e-10    6.8142e-12    4.8479e-07    6.9211e-13
     550    4.3108e-08    7.4433e-10     2.222e-05    2.5091e-14
     600    1.0462e-06    3.6869e-08    0.00053712    2.8358e-11
     650    1.5545e-05    9.9428e-07      0.007892    1.1102e-15
     700    0.00015709    1.6636e-05      0.075039    1.5279e-09
     750     0.0011661    0.00019005       0.39419    4.5887e-08
     800      0.006738     0.0015929       0.70787    6.8612e-12
     850      0.031673      0.010348             1      -0.38394
     900       0.12536      0.054386       0.87501      -0.43902
     950       0.42929       0.23915       0.79689      -0.47807
    1000        1.2998       0.90391       0.79688      -0.49342
    1050        3.5415        3.0015             1      -0.49836
    1100        8.8088        8.9134             1      -0.49952
    1150        20.241        24.024             1       -0.4999
    1200        43.396        59.491       0.38673      -0.49984
    1250         87.53        136.75       0.28907       -0.4998
    1300        167.28        294.36       0.25978      -0.49995
    1350         304.7        597.69       0.24317      -0.49999
    1400        531.75        1152.1       0.21449      -0.49999
    1450        893.03        2119.6       0.25001      -0.49997
    1500        1448.8        3739.6         1e-05          -0.5

Clearly, fsolve appears to get more accurate as ‘T’ gets higher.

The plot:

figure
plot(Tvals, xval)
xlabel('T')
ylabel('x')
grid

looks a bit strange.

[Confidential Confidential]

That's still close to what I'm seeking and I think I can make it work from here. Seriously thank you so much star strider I have been tearing my hair out for three days trying things to get this to work! I can't thank you enough

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.