Trouble finding solution of unknown variable in difficult integration problem

Question

Louis McDermott 2021-1-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/730943-trouble-finding-solution-of-unknown-variable-in-difficult-integration-problem

评论： Louis McDermott 2021-1-31

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2021-1-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/730943-trouble-finding-solution-of-unknown-variable-in-difficult-integration-problem#answer_610218

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Alan Stevens 2021-1-31

编辑：Alan Stevens 2021-1-31

在 MATLAB Online 中打开

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.

Alan Stevens 2021-1-31

在 MATLAB Online 中打开

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 2021-1-31

Yes, this is perfect thank you!

请先登录，再进行评论。

Answer 2

Alex Sha 2021-1-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/730943-trouble-finding-solution-of-unknown-variable-in-difficult-integration-problem#answer_610178

Hi, there are two solutions as below

:

1: h=0.073660616949704

2: h=0.165217787421255

2 个评论
显示无隐藏无

Louis McDermott 2021-1-30

How did you find this?

Alex Sha 2021-1-31

在 MATLAB Online 中打开

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
0.073660616949704
5	0.0711236790825779
0.0690063620895268
5	0.0672422418190603
0.0657674528467493
5	0.0645275378178806
0.0634782731873953

Solution 2:

s	h
0.165217787441994
5	0.156485557183266
0.14921049827769
5	0.143711115798623
0.139599112575786
5	0.136470016301514
0.134031387697928

请先登录，再进行评论。

Trouble finding solution of unknown variable in difficult integration problem

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Trouble finding solution of unknown variable in difficult integration problem

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

2 个评论
显示无隐藏无