LAST EDIT Posted a link to another question with the same problem. Some things should not be done in matlab... maybe someone else will ask this so my conclusion is in the link above.
My coding is hard to read for me and I often make mistakes, so I think a function would help.
EDIT 2 moved old question to the end of it all and improved it
I would like vpasolve to at least automatically add equations or remove them if the unknown is relevant to be calculated.
So let's say I have 4 unknowns: r1, r2, r3, r4 and 4 equations. How can this be done:
syms if r1<=0
r1=0
else
'r1'
same for the other `r`'s;
[if r1<=0
r1=0
else
'r1'
same for the other `r`'s] = vpasolve([if r1<=0
` `
else
'equation1'
same for the other `r`'s], [if r1<=0
r1=0
else
'r1'
same for the other `r's]);
The condition might not always be as easy as >=0...
Also better ways to do it are welcomed as it takes some time to do all this (moved at the end) usingfor j=2:~40. A big time problem could be solved if I also understood why aren't all of my CPU cores used.
Thank you, Dan.
MOVED CONTENT
I would like vpasolve to solve n equations, based on a given n. My n can be 12 or 16, also in the future there will be more than these 2 possibilities. So something like:
function y=custvpas(a,b,c,d,e,f,g,h,n)
if n=12
...
else
...
end
EDIT 1: I have deleted the original vpasolve code as this one is far better and I have added the way it works now.
om(1)=2.76;
om(2)=27.5/56;
om(3)=23.5/56;
om(4)=1;
om(5)=5/1000;
om(6)=50/1000;
om(7)=1/1000;
vb0max=uint8((om(1)+om(3))*10);
cct(1)=om(1)*om(5)/om(6);
c=zeros(15,vb0max);
ct=zeros(4,vb0max);
v=zeros(4,vb0max);
a=zeros(3,vb0max);
v(1,1)=om(6);
ct(1,1)=2*cct(1);
ct(2,1)=om(2)*om(5)/v(1,1);
ct(3,1)=om(3)*om(5)/v(1,1);
cct(4)=cct(1)+ct(2,1)+3/2*ct(3,1);
t = animatedline;
axis([0,vb0max+3,0,7])
for j=1:1
func12
v(1,j+1)=v(1,j)+om(7);
v(2,j+1)=v(2,j)+om(7);
c(14,j)=10^(-14)/c(1,j);
pH(j)=-log10(c(1,j));addpoints(t,v(2,j)*1000,pH(j));drawnow
pause(0.05)
end
for j=2:vb0max+3
v(1,j)=v(1,j-1)+om(7);v(2,j)=v(2,j-1)+om(7);c(15,j)=v(2,j)*om(4)/v(1,j);
if 6*10^(-38)>(c(3,j-1)*v(1,j-1)-0.0001)/v(1,j)*(10^(-14)/(c(1,j-1)*v(1,j-1)-0.0007)/v(1,j))^3
ct(1,j)=(ct(1,j-1)*v(1,j-1)-om(4)*om(7))/v(1,j);
ct(2,j)=ct(2,j-1)*v(1,j-1)/v(1,j);
ct(3,j)=ct(3,j-1)*v(1,j-1)/v(1,j);
v(3,j)=om(7);v(4,j)=0;
func12
c(14,j)=10^(-14)/c(1,j);
else
ct(2,j)=ct(2,j-1)*v(1,j-1)/v(1,j);
func16
if isequal(size(r1), size([0]))
ct(1,j)=(ct(1,j-1)*v(1,j-1)-om(4)*r111)/v(1,j);
ct(3,j)=(ct(3,j-1)*v(1,j-1)-om(4)*r333)/v(1,j);
v(3,j)=r111;v(4,j)=r333;
else
ct(1,j)=(ct(1,j-1)*v(1,j-1)-om(4)*om(7))/v(1,j);
ct(3,j)=ct(3,j-1)*v(1,j-1)/v(1,j);
v(3,j)=om(7);v(4,j)=0;
func12
end
c(14,j)=10^(-14)/c(1,j);
end
pH(j)=-log10(c(1,j));addpoints(t,v(2,j)*1000,pH(j));drawnow
pause(0.05);
end
this is my func12
r1t=ct(1,j);r2t=ct(2,j);r3t=ct(3,j);
syms r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12;
[r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12] = vpasolve([r1t==r1+r5+r6+r8,
r2t==r2+r6+r7,
r3t==r3+r8+r9+r10+r11+r12,
cct(4)==r4+r5+r6+r7+r8+r9+2*r10,
r5==r4*r1*10^1.98,
r6==r4*r1*r2*10^1.08,
r7==r4*r2*10^2.25,
r8==r4*r1*r3*10^2.48,
r9==r4*r3*10^4.04,
r10==r4^2*r3*10^5.38,
r11==1.9*10^(-3)*r3/r1,
r12==2.5*10^(-3)*r3/r1^2], [r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12]);
index = find(r12>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r11>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r10>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r9>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r8>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r7>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r6>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r5>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r4>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r3>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r2>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);index = find(r1>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);
c(1,j)=r1;c(2,j)=r2;c(3,j)=r3;c(4,j)=r4;c(5,j)=r5;c(6,j)=r6;c(7,j)=r7;c(8,j)=r8;c(9,j)=r9;c(10,j)=r10;c(11,j)=r11;c(12,j)=r12;
this is my func16
r2t=ct(2,j);
syms r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 r1t r3t r111 r333;[r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r1t,r3t,r111,r333]=vpasolve([
r1t==r1+r5+r6+r8,
r2t==r2+r6+r7,
r3t==r3+r8+r9+r10+r11+r12,
cct(4)==r4+r5+r6+r7+r8+r9+2*r10,
r5==r4*r1*10^1.98,
r6==r4*r1*r2*10^1.08,
r7==r4*r2*10^2.25,
r8==r4*r1*r3*10^2.48,
r9==r4*r3*10^4.04,
r10==r4^2*r3*10^5.38,
r11==1.9*10^(-3)*r3/r1,
r12==2.5*10^(-3)*r3/r1^2
6*10^(-38)==r3*(10^(-14)/r1)^3,
r1t==(ct(1,j-1)*v(1,j-1)-om(4)*r111)/v(1,j),
r3t==(ct(3,j-1)*v(1,j-1)-om(4)/3*r333)/v(1,j),
om(7)==r111+r333
], [r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r1t,r3t,r111,r333]);
index = find(r10>=0); r10 = r10(index);r1 = r1(index);r2 = r2(index);r3 = r3(index);r4 = r4(index);r5 = r5(index);r6 = r6(index);r7 = r7(index);r8 = r8(index);r9 = r9(index);r11 = r11(index);r12 = r12(index);r111 = r111(index);r333 = r333(index);r1t = r1t(index);r3t = r3t(index);.....r12 = r12(index);r111 = r111(index);r333 = r333(index);r1t = r1t(index);r3t = r3t(index);
