fsolve with two variables in a loop
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi, I thank you in advance.
I try to solve two equation with two variables (x and y) with a changing parameters T, i try to use a fsolve method. I met a problem in the loop. Below is my code.
clear all
% effective normal stress (Pa)
p.sigma=5e7;p.To=273.15+20;p.R=8.314;p.K=6e10;p.Lt=8e-4;p.lambda=0.0625;
p.phio=0.32;p.dsb=2e-6;p.dbulk=2e-5;p.qsb=0.4;p.phic=0.2;p.phio=0.02;
p.qbulk=0.7;p.H=0.577;p.n=1.7;p.m=3;p.An=0.0759;p.At=4.4*p.An;
p.Ea=213e3;p.omega=3.69e-7;p.muo=0.6;p.ao=0.006;p.xo = 0.1813;
p.V0=1e-6;p.V1=0.1e-6;p.gammao=0.02;p.M = 2;p.N = 1;
% a changing parameter T with two variables x and y
T=linspace(273.15,600,100);
x=linspace(0,0.2,100);
y=linspace(0,0.2,100);
% Equation 1: At reference (prestep) velocity of V0
fx=@(x,T) p.Lt.*p.lambda./p.V0*p.An.*exp(-p.Ea./p.R./T).*...
(p.sigma).^p.n./(p.dsb).^p.m.*((p.phic-x)./(x-p.phio)).^(-p.M)-p.H.*(p.qsb-2.*x).^p.N;
% Equation 2: At a velocity (poststep) of V1
fy=@(Y,T) p.Lt.*p.lambda./p.V1*p.An.*exp(-p.Ea./p.R./T).*...
(p.sigma).^p.n./(p.dsb).^p.m.*((p.phic-y)./(y-p.phio)).^(-p.M)-p.H.*(p.qsb-2.*y).^p.N;
musi=@(T) p.muo+p.ao*T./p.To.*log(p.V1./p.V0);
muso=p.muo;
tanso=@(x) p.H.*(p.qsb-2.*x).^p.N;
tansi=@(y) p.H.*(p.qsb-2.*y).^p.N;
aa=@(x,T) p.ao.*T./p.To.*(1+tanso(x).^2)./(1-muso.*tanso(x))./(1-musi(T).*tanso(x));
bb=@(y,T) -tansi(y).*(1+musi(T).^2)./(1-musi(T).*tansi(y))./...
(1-musi(T).*tansi(y)).*(1-(p.V1./p.V0).^(1/(p.M+p.N)))./log(p.V1./p.V0);
% loop
options=odeset('Refine',1,'RelTol',1e-15,'InitialStep',1e-6,'MaxStep',3e6);
for j=1:length(T)
fun=@(x,y)[fx(x,T(j));fy(y,T(j))];
z(:,j)=fsolve(fun,[0.199;0.199],options); %plotted to get initial guess close.
end
% figure
plot(T-273.15,aa(z,T),'k-',T-273.15,aa(z,T),'r-');
2 commentaires
Torsten
le 30 Jan 2023
If the solutions from the symbolic code did not satisfy your constraints, using "fsolve" won't help.
The three solutions you got for x and y from the symbolic approach are the only ones your equations have.
Réponse acceptée
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Stress and Strain dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!