Solution of Non linear Equation

5 ビュー (過去 30 日間)
Monirul Hasan
Monirul Hasan 2019 年 2 月 1 日
回答済み: Monirul Hasan 2019 年 2 月 3 日
Hi, I am using trying to find a sloution for a nonlinear system. With the matlab function 'vpasolve', I am now able to get the solutions but when I am plugging back those solution inside the equation again to check, it seems like they are not the solution for the equations. Is there a way to fix this issue? Code is given below-
clear all, close all, clc
% NonLinear Part rate equation at steady state
syms Ns Nt
assume(Ns > 0);
assume(Nt > 0);
k_rs = 3.2e7 ;
k_ISC = 3.1e7 ;
k_RISC = 5.6e3 ;
k_ST = 2e-10 ;
k_TT = 5e-15 ;
k_NRT = 8.2e2 ;
d = 15e-7 ;
e = 1.6e-19*1e3 ;
J = [0.001:0.5:100];
sol_Ns = zeros(200,3);
sol_Nt = zeros(200,3);
for i = 1:length(J)
eq1 = -(k_rs+k_ISC)*Ns+k_RISC*Nt-k_ST*Ns*Nt+0.25*k_TT*Nt.^2+J(i)/(4*d*e)==0 ;
eq2 = k_ISC*Ns-(k_RISC+k_NRT)*Nt-1.25*k_TT*Nt.^2+(3*J(i))/(4*d*e)==0 ;
sol = vpasolve([eq1, eq2],[Ns,Nt]);
sol_Ns(i,:) = sol.Ns;
sol_Nt(i,:) = sol.Nt;
end

採用された回答

Stephan
Stephan 2019 年 2 月 1 日
編集済み: Stephan 2019 年 2 月 1 日
Hi,
you convert symbolic solutions to double data type. This increases calculation speed, but costs accuracy. Try:
% NonLinear Part rate equation at steady state
syms Ns Nt
assume(Ns > 0);
assume(Nt > 0);
k_rs = 3.2e7 ;
k_ISC = 3.1e7 ;
k_RISC = 5.6e3 ;
k_ST = 2e-10 ;
k_TT = 5e-15 ;
k_NRT = 8.2e2 ;
d = 15e-7 ;
e = 1.6e-19*1e3 ;
J = [0.001:0.5:100];
test = vpa(nan(numel(J),2));
for i = 1:numel(J)
eq1 = -(k_rs+k_ISC)*Ns+k_RISC*Nt-k_ST*Ns*Nt+0.25*k_TT*Nt.^2+J(i)/(4*d*e)==0 ;
eq2 = k_ISC*Ns-(k_RISC+k_NRT)*Nt-1.25*k_TT*Nt.^2+(3*J(i))/(4*d*e)==0 ;
sol(i) = vpasolve([eq1, eq2],[Ns,Nt]);
test(i,1) = -(k_rs+k_ISC)*sol(i).Ns+k_RISC*sol(i).Nt-k_ST*sol(i).Ns*sol(i).Nt+0.25*k_TT*sol(i).Nt.^2+J(i)/(4*d*e);
test(i,2) = k_ISC*sol(i).Ns-(k_RISC+k_NRT)*sol(i).Nt-1.25*k_TT*sol(i).Nt.^2+(3*J(i))/(4*d*e);
end
Best regards
Stephan
  1 件のコメント
Stephan
Stephan 2019 年 2 月 1 日
編集済み: Stephan 2019 年 2 月 1 日
"Hi Stephan, I did'nt understand the new process quite good (not converting symbolic solutions into double data type). Please explain me, as I have just started to learn coding."
There are several data types in Matlab. You are dealing with symbolic variables, which gives exactly results but is very slow. The code you provided in your question contains a conversion form symbolic to the type double. Due to the conversion there is a loss of accuracy, which you see, when you insert the results into the equations. Since you have very big numbers, there are relativly big deviations from zero, which is the expected result, if the solutions are correct. To avoid this i ensured, that the calculation of test keeps symbolic.
"By the way if I want to plot the the solutions as function of J (which is one of my goal), how may I able to do it? Thanks."
To plot the results i changed the code a bit:
% NonLinear Part rate equation at steady state
syms Ns Nt
assume(Ns > 0);
assume(Nt > 0);
k_rs = 3.2e7 ;
k_ISC = 3.1e7 ;
k_RISC = 5.6e3 ;
k_ST = 2e-10 ;
k_TT = 5e-15 ;
k_NRT = 8.2e2 ;
d = 15e-7 ;
e = 1.6e-19*1e3 ;
J = 0.001:0.5:100;
% test = vpa(nan(numel(J),2));
sol_Ns = vpa(nan(numel(J),1));
sol_Nt = vpa(nan(numel(J),1));
for k = 1:numel(J)
eq1 = -(k_rs+k_ISC)*Ns+k_RISC*Nt-k_ST*Ns*Nt+0.25*k_TT*Nt.^2+J(k)/(4*d*e)==0 ;
eq2 = k_ISC*Ns-(k_RISC+k_NRT)*Nt-1.25*k_TT*Nt.^2+(3*J(k))/(4*d*e)==0 ;
[sol_Ns(k), sol_Nt(k)] = vpasolve([eq1, eq2],[Ns,Nt]);
% test(k,1) = -(k_rs+k_ISC)*sol_Ns(k)+k_RISC*sol_Nt(k)-k_ST*sol_Ns(k)*sol_Nt(k)+0.25*k_TT*sol_Nt(k).^2+J(k)/(4*d*e);
% test(k,2) = k_ISC*sol_Ns(k)-(k_RISC+k_NRT)*sol_Nt(k)-1.25*k_TT*sol_Nt(k).^2+(3*J(k))/(4*d*e);
end
yyaxis left
plot(J, sol_Ns)
yyaxis right
plot(J, sol_Nt)
Now the results are plotted as you wanted to do. I commented the part of test out, since we saw that the solutions are correct one time. If you need it just uncomment the 3 lines regarding test.
I also changed the counter variable from 'i' to 'k', since 'i' is used for complex numbers. It is good practice to avoid 'i' as counter variable.

サインインしてコメントする。

その他の回答 (2 件)

Monirul Hasan
Monirul Hasan 2019 年 2 月 1 日
Hi Stephan, I did'nt understand the new process quite good (not converting symbolic solutions into double data type). Please explain me, as I have just started to learn coding. By the way if I want to plot the the solutions as function of J (which is one of my goal), how may I able to do it? Thanks.

Monirul Hasan
Monirul Hasan 2019 年 2 月 3 日
Thanks a lot Stephan.

カテゴリ

Help Center および File ExchangeRobust Control Toolbox についてさらに検索

タグ

製品


リリース

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by