Trouble finding solution of unknown variable in difficult integration problem

I am trying to solve for 'h' using the first equation above. My script produces a result, however it is unexpected. The solution of 'h' should increase with 's', yet this is not the case. I have attached the script below, any help would be much appreciated.

clc; clear all; close all;

%%%%%%%%%% Physical Parameters %%%%%%%%%%

r=0.1;  
b=0.15;
W=60;

Kc=0.99*1000;    
Kphi=1528.4*1000;   
n=1.1;   
c=1.04*1000;    
phi=28; 

A1=(Kc/b+Kphi);

k=0.6;

Beta=0.0872665;

kx=0.043*Beta+0.036;

%%%%%%%%%% Defining equations %%%%%%%%%%

s=1
    
syms h 
syms Pheta

PhetaF=acos(1-h/r);
PhetaR=acos(1-k*h/r);

jx=r*(PhetaF-Pheta-(1-s)*(sin(PhetaF)-sin(Pheta)));

Sig=A1*r^n*(cos(Pheta)-cos(PhetaF));

Taux=(c+Sig*tand(phi))*(1-exp(-jx/kx));

%%%%%%%%%% Solving for h %%%%%%%%%%

fun=Taux*sin(Pheta)+Sig*cos(Pheta);

eqn=W==r*b*vpaintegral(fun,Pheta,[PhetaR,PhetaF]);

sol=vpasolve(eqn,h,.1)

Answer 1

[Alan Stevens]

I used a different approach - see below. I could only see an increase in h, with increasing s, as a step change around s = 4.5.

s = 1:0.1:10;
h0 = 0.1;   % Use h0 = 0.15 to get a different set of results.
for i = 1:numel(s)
    h(i) = fzero(@(h) Zfun(h,s(i)), h0);
end

plot(s,h,'o'),grid
xlabel('s'),ylabel('h')
legend(['initial guess h0 = ' num2str(h0)])
disp(h)

function Z = Zfun(h,s)

        r=0.1;  
        n=1.1;  
        rn = r^n;
        b=0.15;
        W=60;
        Kc=0.99*1000;    
        Kphi=1528.4*1000;
        k = 0.6;
        c=1.04*1000;    
        tanphi=tand(28); 
        Beta=0.0872665;
        kx=0.043*Beta+0.036;
        sigA = (Kc/b+Kphi)*rn;

        thetar = acos(1-k*h/r);
        thetaf = acos(1-h/r);
        
        sig = @(theta) sigA*(cos(theta) - cos(thetaf));
        jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
        taux = @(theta) (c + sig(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
        kern = @(theta) taux(theta).*sin(theta) + sig(theta).*cos(theta);
        Z = r*b*integral(kern,thetar,thetaf) - W;
        
end

This results in

[Louis McDermott]

Would it be possible to use these values of h to calculate new values of thetaF and thetaR and use them to calculate:

Fx=r*b*int(taux1(theta).*cos(theta)-sig1(theta).*cos(theta),thetaR,thetaF)

I know h is now an array, and so it is impossible to integrate over array boundaries, is there a way around this?

[Alan Stevens]

Do you mean something like this:

s = 1:0.1:10;
h0 = 0.1;   % Use h0 = 0.15 to get a different set of results.
for i = 1:numel(s)
    h(i) = fzero(@(h) Zfun(h,s(i)), h0);
    Fx(i) = integralfunction(h(i),s(i));
end

plot(s,h,'o'),grid
xlabel('s'),ylabel('h')
legend(['initial guess h0 = ' num2str(h0)])
%disp(h)
figure
subplot(2,1,1)
plot(h,Fx,'o'),grid
xlabel('h'),ylabel('Fx')
subplot(2,1,2)
plot(s,Fx,'o'),grid
xlabel('s'),ylabel('Fx')


function Z = Zfun(h,s)
        W=60;
        Z = integralfunction(h,s) - W;
end

function Irb = integralfunction(h,s)

       r=0.1;  
        n=1.1;  
        rn = r^n;
        b=0.15;
       
        Kc=0.99*1000;    
        Kphi=1528.4*1000;
        k = 0.6;
        c=1.04*1000;    
        tanphi=tand(28); 
        Beta=0.0872665;
        kx=0.043*Beta+0.036;
        sigA = (Kc/b+Kphi)*rn;

        thetar = acos(1-k*h/r);
        thetaf = acos(1-h/r);
        
        sig = @(theta) sigA*(cos(theta) - cos(thetaf));
        jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
        taux = @(theta) (c + sig(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
        kern = @(theta) taux(theta).*sin(theta) + sig(theta).*cos(theta);
        Irb = r*b*integral(kern,thetar,thetaf);

end

Incidentally, did you mean to put a negative sign between the two terms in the integral? I've kept the positive sign in the above to match the original.

[Louis McDermott]

This function, Fx, is a new function which uses the now known values of thetar and thetaf using the calculated h, hence the minus should be there, also the cos and sin have switched.

[Alan Stevens]

Ok. Assuming you want the original values of h, then you can do the following:

s = 1:0.1:10;
h0 = 0.1;   % Use h0 = 0.15 to get a different set of results.
for i = 1:numel(s)
    h(i) = fzero(@(h) Zfun(h,s(i)), h0);
    Fx(i) = integralfunction2(h(i),s(i));
end

plot(s,h,'o'),grid
xlabel('s'),ylabel('h')
legend(['initial guess h0 = ' num2str(h0)])
%disp(h)
figure
subplot(2,1,1)
plot(h,Fx,'o'),grid
xlabel('h'),ylabel('Fx')
subplot(2,1,2)
plot(s,Fx,'o'),grid
xlabel('s'),ylabel('Fx')


function Z = Zfun(h,s)
        W=60;
        Z = integralfunction(h,s) - W;
end

function Irb = integralfunction(h,s)

       r=0.1;  
        n=1.1;  
        rn = r^n;
        b=0.15;
       
        Kc=0.99*1000;    
        Kphi=1528.4*1000;
        k = 0.6;
        c=1.04*1000;    
        tanphi=tand(28); 
        Beta=0.0872665;
        kx=0.043*Beta+0.036;
        sigA = (Kc/b+Kphi)*rn;

        thetar = acos(1-k*h/r);
        thetaf = acos(1-h/r);
        
        sig = @(theta) sigA*(cos(theta) - cos(thetaf));
        jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
        taux = @(theta) (c + sig(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
        kern = @(theta) taux(theta).*sin(theta) + sig(theta).*cos(theta);
        Irb = r*b*integral(kern,thetar,thetaf);

end


function Irb2 = integralfunction2(h,s)
       r=0.1;  
        n=1.1;  
        rn = r^n;
        b=0.15;
       
        Kc=0.99*1000;    
        Kphi=1528.4*1000;
        k = 0.6;
        c=1.04*1000;    
        tanphi=tand(28); 
        Beta=0.0872665;
        kx=0.043*Beta+0.036;
        sigA = (Kc/b+Kphi)*rn;

        thetar = acos(1-k*h/r);
        thetaf = acos(1-h/r);
        
        sig1 = @(theta) sigA*(cos(theta) - cos(thetaf));
        jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
        taux1 = @(theta) (c + sig1(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
        kern1 = @(theta) taux1(theta).*cos(theta)-sig1(theta).*cos(theta);

        Irb2 = r*b*integral(kern1,thetar,thetaf);
end

[Louis McDermott]

Yes, this is perfect thank you!

Answer 2

[Alex Sha]

Hi, there are two solutions as below:

1: h=0.073660616949704

2: h=0.165217787421255

[Louis McDermott]

How did you find this?

[Alex Sha]

If give s form 1 to 4, there will be two set of solutions, however, both of them, 'h' will decrease with 's':

solution 1:

s	h
1	0.073660616949704
1.5	0.0711236790825779
2	0.0690063620895268
2.5	0.0672422418190603
3	0.0657674528467493
3.5	0.0645275378178806
4	0.0634782731873953

Solution 2:

s	h
1	0.165217787441994
1.5	0.156485557183266
2	0.14921049827769
2.5	0.143711115798623
3	0.139599112575786
3.5	0.136470016301514
4	0.134031387697928