fsolve with two variables in a loop
Show older comments
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-');
Accepted Answer
More Answers (0)
Categories
Find more on Programming 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!