Solution of Non linear Equation

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

Answer 1

[Stephan]

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

[Stephan]

"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.

Answer 2

[Monirul Hasan]

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.

Answer 3

[Monirul Hasan]

Thanks a lot Stephan.