Solving nonlinear equation using fsolve.

2015 11 30
1 个回答
更新时间：2015 12 1
4 次查看（30 天）

0 个投票

I have nonlinear coupled equation. I have T varied from 300 to 1300, each T give coupled-equation. How to solve each equation using fsolve? I know to solve coupled equation if only there is one coupled-equation, but do I have to solve one by one? 1000 equation total. Thank you

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

回答（1 个）

Torsten
Torsten 2015-11-30

0 个投票

Use a for-loop in which you call fsolve (or fzero) 1000 times with T varying from 300 to 1300.
Best wishes
Torsten.

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

Yes, I have tried using for-loop, but the function must be saved in another file so i can call it. but, do i have to save 1000 function/equation? as you can see in:
http://www.mathworks.com/help/optim/ug/fsolve.html?searchHighlight=fsolve#butbm05-1
and
http://www.dam.brown.edu/people/sahn/html/fsolveDemo.html
Can I add fsolve directly in my loop? pls help me with the coding. thank you.
Torsten
Torsten 2015-12-1
编辑：Torsten 2015-12-1
If I understand you correctly, the equation you are trying to solve vary with T.
Thus my suggestion is to call fsolve somehow like
x0 = 1;
for i = 1:1001
T = 300+(i-1);
sol(i) = fsolve(@(x)fun(x,T),x0);
x0 = sol(i);
end
This way, you get a solution array sol(i) depending on the temperature T(i).
Best wishes
Torsten.
kusniar deny permana
编辑：kusniar deny permana 2015-12-1
Here is my equation and code:
syms x
%R,Regnault constant (gas constant) = 8.3144598 J/(K.mol)
R=8.3144598
%temperature
T=300
Z=1
while T<=1300
%Ga1, Gibbs standard for Mg alpha(hcp)
if T>=298 & T<=923
Ga1(Z)= (-8367.34)+(143.675547.*T)-(26.1849782.*T.*log(T))+(0.4858*(10.^-3).*(T.^2))-(1.393669*10.^-6.*T.^3)+(78950.*T.^-1)
else T>923 & T<=1300
Ga1(Z)=(-14130.185)+(204.716215.*T)-(34.3088.*T.*log(T))+(1038.192.*10.^25.*T.^-9)
end
%Gb1, Gibbs standard for Li alpha(hcp)
if T>=200 & T<=453
Gb1(Z)=(-10737.817)+(219.637482.*T)-(38.940488.*T.*log(T))+(35.466931.*10.^-3.*T.^2)-(19.869816.*10.^-6.*T.^3)+(159994.*T.^-1)
elseif T>453 & T<=500
Gb1(Z)=(-559733.123)+(10549.879893.*T)-(1702.8886493.*T.*log(T))+(2258.3294448.*10.^-3.*T.^2)-(571.066077.*10.^-6.*T.^3)+(33885874.*T.^-1)
else T>500 & T<=1300
Gb1(Z)=(-9216.994)+(181.278285.*T)-(31.2283718.*T.*log(T))+(2.633221.*10.^-3.*T.^2)-(0.438058.*10.^-6.*T.^3)-(102387.*T.^-1)
end
%Ga2, Gibbs standard for Mg beta(bcc)
if T>=298 & T<=923
Ga2(Z)=(-5267.34)+(141.575547.*T)-(26.1849782.*T.*log(T))+(0.4858.*10.^-3.*T.^2)-(1.393669.*10.^-6.*T.^3)+(78950.*T.^-1)
else T>923 & T<=1300
Ga2(Z)=(-11030.185)+(202.616215.*T)-(34.3088.*T.*log(T))+(1038.192*10.^25.*T.^-9)
end
%Gb2, Gibbs standard for Li beta(bcc)
if T>=200 & T<=453
Gb2(Z)=(-10583.817)+(217.637482.*T)-(38.940488.*T.*log(T))+(35.466931.*10.^-3.*T.^2)-(19.869816.*10.^-6.*T.^3)+(159994.*T.^-1)
elseif T>453 & T<=500
Gb2(Z)=(-559579.123)+(10547.879893.*T)-(1702.8886493.*T.*log(T))+(2258.329444.*10.^-3.*T.^2)-(571.066077.*10.^-6.*T.^3)+(33885874.*T.^-1)
else T>500 & T<=1300
Gb2(Z)=(-9062.994)+(179.278285.*T)-(31.2283718.*T.*log(T))+(2.633221.*10.^-3.*T.^2)-(0.438058.*10^-6.*T.^3)-(102387.*T.^-1)
end
%Xb=mol fraction Li, Xa=mol fraction Mg, x defined as Xb so Xa=(1-x)
Gm1= (1-x).*Ga1+x.*Gb1+R.*T.*((1-x).*log(1-x)+x.*log(x))+(1-x).*x.*((-6856)+4000.*(1-2.*x)+4000.*(1-2.*x).^2)
Gm2= (1-x).*Ga2+x.*Gb2+R.*T.*((1-x).*log(1-x)+x.*log(x))+(1-x).*x.*((-18335)+8.49.*T+3481.*(1-2.*x)+(2658-0.114.*T).*(1-2.*x).^2)
%assumption Gm1=Gm2 to find roots x0, intersection two phase alpha-beta
%G1 is difference of Gibbs Li and Mg
Gm=Gm2-Gm1
%roots of Gm, x0
for i=1:length(Gm)
Gma=Gm(i)
Gma==0
Gmb=solve(Gma,x)
Gmd(i)=Gmb(2)
end
F1=diff(Gm1)
F2=diff(Gm2)
F = [F1-F2; %in x1 and x2
F1.*(x(2)-x(1))-(Gm2(x(2))-Gm1(x(1)))];
% loop for next T value
T=T+100
Z=Z+1
end
The equation I want to solve is coupled equation "F", which is roots x1 and x2. I found x0 as Gmd. So i have a coupled equation and x0 as T change.
I need to solve coupled-equation as T change (on while-loop). So it would be Ti,x1i,x2i.......Tn,x1n,x2n. Thank you Torsten.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by