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);

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.

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.

Star Strider
on 19 Nov 2020

My pleasure!

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

.

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

Start Hunting!
## 0 Comments

Sign in to comment.