system of non linear equations

Hello

I am dealing with a system of two non-linear equations and I tried to solve it by 'fsolve', but the solver stopped and could not fine the answers.

Would you please help me?

I attached the files

Answer 1

[Torsten]

%solving stsyem of nonlinear equation 3.31 and 3.32 Woods paper with fsolve
% eta_l= x1 , eta_u= x2
r=2;
v=20;
F=13/9;

fun= @(x)rooted(x,r,v,F);
x0 = [0 0];
x = fsolve(fun,x0)

No solution found.

fsolve stopped because the problem appears regular as measured by the gradient,
but the vector of function values is not near zero as measured by the
value of the function tolerance.

x = 1×2

0 0

function f = rooted (x,r,v,F)
  %UNTITLED3 Summary of this function goes here
  %   Detailed explanation goes here
  % F<Fd Fd= 2.6 I put F=1.5
  f(1)= (r*v-1)/(v*(r-1)) * ((r*v*x(1)^2-x(2)^2)/(r*v-1))^(3/2) - ...
        (v-1)/(v*(r-1)) * ((v*x(1)^2-x(2)^2)/(v-1))^(3/2) - 9 ;
  f(2)= x(2)^3/(r*v) + (v-1)/(v*(r-1)) * ((v*x(1)^2-x(2)^2)/(v-1))^(3/2) - ... 
        (r*v-1)/(r*v*(r-1)) * ((r*v*x(1)^2-x(2)^2)/(r*v-1))^(3/2) - 9*F;
end

[sepideh]

yes I received the same answer.

Do you have any suggestion how to solve thses equations?

[Torsten]

Do you have any suggestion how to solve thses equations?

Yes, by using better initial guesses for the solution than [0 0].

Your implementation of the equations had several errors - so better use the corrected code from above when trying further.

[sepideh]

sure.

Thank you.

I will try with corrected code

[sepideh]

I checked with corrected code and new initial guess .it worked now, Thank you

[sepideh]

Dear Torsten

I have the similar problem. this time the non linear system is a bit different.

would you please check this too?

again the solver can not find the roots

I put different initial guess but still no results

thank you in advance

% upper layer detacheded 

% plane geometry 

clc
clear all
close all

%solving stsyem of nonlinear equation 3.38 and 3.39 Woods paper with fsolve
% eta_l= x1 , eta_u= x2

%%%% the condition r*v *xil^2>xiu^2> xil^2 must be held  %%%


r=2;
v=0.1;
 Fu=  1/( ( (1-v)/(v*(r-1)))^(1/2)-1)   %eq 3.37 is used here as indicator 

Fu = 0.5000

% gam > Fu



F=1.44;
fun= @(x)rooted(x,r,v,F);
x0 = [0 5];
x = fsolve(fun,x0)

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

x = 1×2

0 1.4104

fun(x)

ans = 1×2

-7.2060 0.9307

xil= x(1)

xil = 0

xiu=x(2)

xiu = 1.4104

xi_star=((xiu^2-v*xil^2)/(1-v))^(1/2) %here is the point that lower layer goes to zero, or upper layer starts

xi_star = 1.4867

xi_plus= ((xiu^2-r*v*(xil^2))/(1-r*v))^(1/2)%here is the point that upper layer goes to zero.

xi_plus = 1.5769

% 0<xi<xi+ is the region that only lower layer is present

xi = linspace(0,xi_plus,1000);
fll= @(xi) ((xiu.^2-xi.^2)./r*v*6);
plot(xi,fll(xi),'r')
hold on 
xii = linspace(xi_plus,xi_star,1000);

fl=@(xii) ((xiu^2-(v*xil^2))/(6*v*(r-1)))+(((v-1)*xii.^2)/(6*v*(r-1)));
fu=@(xii) (((r*v*xil^2)-xiu^2)/(6*v*(r-1)))+(((1-r*v)*xii.^2)/(6*v*(r-1)))+((xiu^2-(v*xil^2))/(6*v*(r-1)))+(((v-1)*xii.^2)/(6*v*(r-1)));


plot (xii,fl(xii),'r')
hold on 
plot (xii,fu(xii),'b')
hold on
xiii=linspace(xi_star,xil,1000);
fuu=@(xiii) ((xil^2-xiii.^2)./(6));
plot (xiii,fuu(xiii),'b')
hold off

B =area(xi, fll(xi));
B.FaceColor = [1 0 0]; % red
B.LineWidth = 1;
hold on
C =area(xii, fu(xii));
C.FaceColor = [0 0 1];
C.EdgeColor = [0 0 1];
C.LineWidth = 1;
hold on
d =area(xii, fl(xii));
d.FaceColor = [1 0 0]; % red
d.EdgeColor = [1 0 0];
d.LineWidth = 1;
hold on
A =area(xiii, fuu(xiii));
A.FaceColor = [0 0 1];
A.EdgeColor = [0 0 1];
A.LineWidth = 1;

xlabel('ξ')
ylabel('f(ξ)')


legend([C B],'Upper Layer','Lower Layer')

function f = rooted (x,r,v,F)
%eq 3.38 and 3.39
  
  % F>Fu Fu= 0.500 I put F=1.5
  f(1)= (1-r*v)/(v*(r-1)) * ((x(2)^2-r*v*x(1)^2)/(1-r*v))^(3/2) - ...
        (1-v)/(v*(r-1)) * ((x(2)^2-v*x(1)^2)/(1-v))^(3/2) + x(1)^3 - 9 ;
  f(2)=  (1-v)/(v*(r-1)) * ((x(2)^2-v*x(1)^2)/(1-v))^(3/2) - ... 
        (1-r*v)/(r*v*(r-1)) * ((x(2)^2-r*v*x(1)^2)/(1-r*v))^(3/2) - 9*F;
end

[Torsten]

I don't know what the reason for non-convergence is. The equations seem to be correct - and concerning parameters r, v and F and initial values for x - I don't know if they make sense.

[sepideh]

I chose r ,v and initial guess based on the paper and it should produces the graph I attached.

lower layer (red one) remained attached to the source, while upper layer (blue one) must over run and detached from the source.