MATLAB Answers

How to find the upper limit of a complex integral

2 views (last 30 days)
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);

  0 Comments

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 19 Nov 2020
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.

  4 Comments

Show 1 older comment
Star Strider
Star Strider on 19 Nov 2020
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
Confidential Confidential on 19 Nov 2020
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
Star Strider on 19 Nov 2020
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by