Im trying to run fsolve inside a for loop but it seem that my second iteration inside the for loop is not able to read the initial trial point
1 view (last 30 days)
Show older comments
Code below:
clear;
clc;
%%
Disp = [1 10 4 2 3;
2 8 3 2 2;
3 5 4 1.5 3];
options = optimset('TolFun',1e-12,'MaxFunEvals',1e9,'Maxiter',1e9);
options.StepTolerance =1e-12;
%%
N = size(Disp,1);
O = [0,0];
tau = 0.3;
Bdispy = 0;
Cdispy = 0;
GG = zeros(2,N);
exitflags = zeros(1,N);
for i = 1:N
Peakdispx(i) = Disp(i,2);
Peakdispy(i) = Disp(i,3);
Adispy(i) = Peakdispy(i);
mslope(i) = (Disp(i,5)-O(2))/(Disp(i,4)-O(1));
Cdispx(i) = Peakdispx(i)-((Peakdispy(i)-Cdispy)/(mslope(i)));
Bdispx(i) = Cdispx(i)*-1;
Adispx(i) = (Adispy(i)-Bdispy)/(mslope(i))+Bdispx(i);
slopeACdisp(i) = ((Adispy(i)-Cdispy)/(Adispx(i)-Cdispx(i)));
Afull(i) = (0.5*abs(Adispx(i)*Bdispy+Bdispx(i)*Cdispy+Cdispx(i)*Peakdispy(i)+Peakdispx(i)*Adispy(i)-Adispy(i)*Bdispx(i)-Bdispy*Cdispx(i)-Cdispy*Peakdispx(i)-Peakdispy(i)*Adispx(i)));
F = @(GG)[slopeACdisp(i)-((GG(2)-Cdispy)/(GG(1)-Cdispx(i)));
-tau+(0.5*abs((GG(1)*Bdispy+Bdispx(i)*Cdispy+Cdispx(i)*Peakdispy(i)+Peakdispx(i)*GG(2)-GG(2)*Bdispx(i)-Bdispy*Cdispx(i)-Cdispy*Peakdispx(i)-Peakdispy(i)*GG(1))))/Afull]
gg0 = [tau*Disp(i,3);(1-tau)*Disp(i,2)];
GG = fsolve(F,gg0);
%% pinch in x and y
disppinchy(i) = GG(2)/Disp(i,3);
U(i) = ((GG(2)-Cdispy)/mslope(i))+Cdispx(i);
L(i) = ((GG(2)-Bdispy)/mslope(i))+Bdispx(i);
disppinchx(i) = (GG(2)-L(i))/(U(i)-L(i));
end
for i = 1:N
Gdispx(i) = GG(1);
Gdispy(i) = GG(2);
Apin(i) =0.5*abs(Gdispx(i)*Bdispy+Bdispx(i)*Cdispy+Cdispx(i)*Peakdispy(i)+Peakdispx(i)*Gdispy(i)-Gdispy(i)*Bdispx(i)-Bdispy*Cdispx(i)-Cdispy*Peakdispx(i)-Peakdispy(i)*Gdispx(i));
ratio(i) = Apin(i)/Afull(i)
end
0 Comments
Answers (1)
Matt J
on 25 Nov 2021
Edited: Matt J
on 25 Nov 2021
Should the division be by Afull(i)?
F = @(GG)[slopeACdisp(i)-((GG(2)-Cdispy)/(GG(1)-Cdispx(i)));
-tau+(0.5*abs((GG(1)*Bdispy+Bdispx(i)*Cdispy+Cdispx(i)*Peakdispy(i)+Peakdispx(i)*GG(2)-GG(2)*Bdispx(i)-Bdispy*Cdispx(i)-Cdispy*Peakdispx(i)-Peakdispy(i)*GG(1))))/Afull(i)]
2 Comments
Matt J
on 25 Nov 2021
Since it's a 2D function, you can do a grid search for a more accurate initial guess.
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!