Why my code find a second intersection of curve with line and not the first?
2 ビュー (過去 30 日間)
古いコメントを表示
I am trying to find a intersection of very complex curve's equation with line parallel with X-axis, but Matlab gives me a result which is suitable for the second intersection. I need first one.
Plot code for better understand is also included. Magenta line should reach and stop at first intersection (pH=6.4 not the second one at pH=13,38). What am I doing wrong?
Oh yes, and I need a exact value of X-coordinate of the first intersection and also a graphical solution as it can be seen.
Here is a MWE:
clc;
close all;
%There are input constants
cP=0.05;
cN=0.05;
cMg=0.005;
logb0=2.88;
logb1=1.17;
logb2=4.85;
logb3=2.58;
logKs1=-12.6;
pKa1=-2.143;
pKa2=-7.205;
pKa3=-12.34;
pKb1NH3=-4.753355133;
%there is a range of x-axis (0-14)
pH=0:0.01:14;
%There are calculations of cH2PO4 cHPO4 cPO4, they needs to be calculated for final curve equation (they are variables of final curve equation and depend on pH value)
Ka1=power(10,pKa1);
Ka2=power(10,pKa2);
Ka3=power(10,pKa3);
cH2PO4=cP.*Ka1.*(power(10,-pH)).^(2)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cHPO4=cP.*Ka1.*Ka2.*power(10,-pH)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cPO4=cP.*Ka1.*Ka2.*Ka3./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
%There is calculations of cNH4, it have to be calculated for final curve equation (it is variable of final curve equation and depends on pH value also)
pOH=14-pH;
Kb1=power(10,pKb1NH3);
cNH4=cN.*Kb1./(power(10,-pOH)+Kb1);
%there are some partial calculations for the final curve equation
b0=power(10,logb0);
b1=power(10,logb1);
b2=power(10,logb2);
b3=power(10,logb3);
Ks1=power(10,logKs1);
cOH=power(10,pH-14);
%this is the final curve equation
logcNP=Ks1./(cPO4.*cNH4).*(1+b1.*cH2PO4+b0.*cHPO4+b2.*cPO4+b3*cOH);
%I am trying to calculate a intersection of the curve above with y value of 0.005 (cMg at the start of the code).
crossectionY=cMg
[crossectionY, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))
%This is plot of the final curve equation (black) and calculated intersection with magenta line
semilogy(pH,logcNP,'-', 'color', 'k', 'MarkerFaceColor','w', 'MarkerEdgeColor', 'm');
box off;
hold on;
plot([0 crossectionX],[cMg cMg],'-', 'color', 'm', 'LineWidth',0.5);
box off;
hold on;
axis([0 14 0.000000000000000001 1000000000000000000]);
xlabel('pH (-)');
ylabel('c {(mol/dm^{3})}');
axis square;
grid on;
axis on;
box off;
hold off;
Than you for answers. :-)
0 件のコメント
採用された回答
Star Strider
2023 年 5 月 14 日
See if the lines oif copde I added do what you want —
%There are input constants
cP=0.05;
cN=0.05;
cMg=0.005;
logb0=2.88;
logb1=1.17;
logb2=4.85;
logb3=2.58;
logKs1=-12.6;
pKa1=-2.143;
pKa2=-7.205;
pKa3=-12.34;
pKb1NH3=-4.753355133;
%there is a range of x-axis (0-14)
pH=0:0.01:14;
%There are calculations of cH2PO4 cHPO4 cPO4, they needs to be calculated for final curve equation (they are variables of final curve equation and depend on pH value)
Ka1=power(10,pKa1);
Ka2=power(10,pKa2);
Ka3=power(10,pKa3);
cH2PO4=cP.*Ka1.*(power(10,-pH)).^(2)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cHPO4=cP.*Ka1.*Ka2.*power(10,-pH)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cPO4=cP.*Ka1.*Ka2.*Ka3./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
%There is calculations of cNH4, it have to be calculated for final curve equation (it is variable of final curve equation and depends on pH value also)
pOH=14-pH;
Kb1=power(10,pKb1NH3);
cNH4=cN.*Kb1./(power(10,-pOH)+Kb1);
%there are some partial calculations for the final curve equation
b0=power(10,logb0);
b1=power(10,logb1);
b2=power(10,logb2);
b3=power(10,logb3);
Ks1=power(10,logKs1);
cOH=power(10,pH-14);
%this is the final curve equation
logcNP=Ks1./(cPO4.*cNH4).*(1+b1.*cH2PO4+b0.*cHPO4+b2.*cPO4+b3*cOH);
%I am trying to calculate a crosssection of the curve above with y value of 0.005 (cMg at the start of the code).
crossectionY=cMg
[crossectionY, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))
xcidx = find(diff(sign(logcNP - 0.005))); % ADDED: Approximate Crossing In dices For 'logcNP' At 0.005
for k = 1:numel(xcidx)
idxrng = max(1,xcidx(k)-1) : min(numel(pH),xcidx(k)+1); % ADDED: Index Range For Interpolation
pHx(k) = interp1(logcNP(idxrng),pH(idxrng), 0.005); % ADDED: Values Of 'pH' For 'logcNP = 0.005'
end
pHx
%This is plot of the final curve equation (black) and calculated crossection with magenta line
semilogy(pH,logcNP,'-', 'color', 'k', 'MarkerFaceColor','w', 'MarkerEdgeColor', 'm');
box off;
hold on;
plot([0 crossectionX],[cMg cMg],'-', 'color', 'm', 'LineWidth',0.5);
plot(pHx, ones(size(pHx))*0.005, 'sc', 'MarkerSize', 10) % ADDED: Plot Cyan Squares At Intersection Points
box off;
hold on;
axis([0 14 0.000000000000000001 1000000000000000000]);
xlabel('pH (-)');
ylabel('c {(mol/dm^{3})}');
axis square;
grid on;
axis on;
box off;
hold off;
.
2 件のコメント
Star Strider
2023 年 5 月 15 日
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Lighting, Transparency, and Shading についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!