Plot multiple graph after solving

I have following two equations -

    f1 = -2*g*t.^3/((1 - y).^(1-2*g)) + 2*g*(3 - t.^2)^2./((2 - t).*y.^(1-2*g)) + (8*g*t)./(y.^(1-2*g)) - 2*g*(3 - 2*t).^2./((2 - t).*(1 - y).^(1-2*g)); 
    f2 = t - (((y.^(2*g))*(2/3) - (((1-y).^(2*g)).*((3 - 2*t)./(3*(2-t)))))./((y.^(2*g)).*(((3-t.^2))./(3*(2-t))) - ((1-y).^(2*g)).*(t/3)));

For every g in the range (0,0.114), solving the two equations above give two solution - y with t=1 and y with a t between 0.5 and 1. (it also give solution with negative y and t less than 0.5, I am not interested in them). To give a sense, for g=0.11 I get following solution - t=1 ^ y=0.85536 t=0.9873 ^ y=0.85530.

I want Matlab to plot two graph (on same plot) with g on x axis and y values in y axis. One graph should have y values for t=1 and one graph should have y values for t in range 0.5 to 1 (below 1).

Since I have done numerical simulation with random values, I know graph should look the following (attached rough sketch)

I don't know how to ask Matlab to pick the specific ans and then plot them. Shall I create a loop?

5 Comments

Can you explain more about how you are solving those equations and, in particular, how the solutions are stored?
I just simply used syms command to simultaneously solve these equation. The main problem is, as you mentioned, that I don't know how to ask Matlab to store the solution for plotting and secondly how to separate the solution (one corresponding to t=1 and other corresponding to t<1) for the two graphs that I want
I have following two equations
You have functions, not equations. Or do you mean f1 = 0 and f2 = 0 ?
Yes, both are equal to zero, sorry for the confusion
@Torsten Please guide me how to proceed. I have just given the two equation name for Matlab to work with. Both f1 and f2 are equal to zero. I can solve them for one particular value of g but don't know how to store the obtained values and to create two graphs into same plot using that data

Sign in to comment.

 Accepted Answer

syms y t g
f1 = -2*g*t.^3/((1 - y).^(1-2*g)) + 2*g*(3 - t.^2)^2./((2 - t).*y.^(1-2*g)) + (8*g*t)./(y.^(1-2*g)) - 2*g*(3 - 2*t).^2./((2 - t).*(1 - y).^(1-2*g));
[N1,D1] = numden(f1);
f2 = t - (((y.^(2*g))*(2/3) - (((1-y).^(2*g)).*((3 - 2*t)./(3*(2-t)))))./((y.^(2*g)).*(((3-t.^2))./(3*(2-t))) - ((1-y).^(2*g)).*(t/3)));
[N2,D2] = numden(f2);
tsol = solve(N2==0,t,'MaxDegree',3);
N1_subs_tsol1 = subs(N1,t,tsol(1));
N1_subs_tsol2 = subs(N1,t,tsol(2));
N1_subs_tsol3 = subs(N1,t,tsol(3));
G = 0.1:0.001:0.110;
soly1 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G)==0),G);
solt1 = arrayfun(@(soly1,G)subs(tsol(1),[y,g],[soly1,G]),soly1,G);
soly2 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G)==0),G);
solt2 = arrayfun(@(soly2,G)subs(tsol(2),[y,g],[soly2,G]),soly2,G);
soly3 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G)==0),G);
solt3 = arrayfun(@(soly3,G)subs(tsol(3),[y,g],[soly3,G]),soly3,G);
hold on
p1 = plot(G,soly1,'r');
%plot(G,solt1,'r')
p2 = plot(G,soly2,'b');
%plot(G,solt2,'b')
p3 = plot(G,1./(1+(1/3).^(1./(1-2*G))),'g');
legend([p1,p2,p3],["slice 1","slice 2","slice 3"])
xlabel("g")
hold off
grid on
% Check the solutions
err11 = arrayfun(@(solt1,soly1,G)double(subs(f1,[t y g],[solt1 soly1 G])),solt1,soly1,G)
err11 = 1×11
1.0e-32 * 0 -0.0000 -0.0000 0.0000 0.0000 0.0000 0 -0.0000 0.0000 0.0000 0.1255
err12 = arrayfun(@(solt1,soly1,G)double(subs(f2,[t y g],[solt1 soly1 G])),solt1,soly1,G)
err12 = 1×11
1.0e-39 * 0 0 0 0 0.1837 0 0 0 0 0.1837 0
err21 = arrayfun(@(solt2,soly2,G)double(subs(f1,[t y g],[solt2 soly2 G])),solt2,soly2,G)
err21 = 1×11
1.0e-33 * 0.0000 0.0007 -0.0000 0.0000 0.9884 -0.0000 0.0001 0.0000 -0.3493 0.0000 -0.0000
err22 = arrayfun(@(solt2,soly2,G)double(subs(f2,[t y g],[solt2 soly2 G])),solt2,soly2,G)
err22 = 1×11
1.0e-39 * 0.3673 -0.5510 0.5510 -0.1837 -0.1837 0 -0.5510 0.1837 0 -0.5510 -0.9184
err31 = arrayfun(@(solt3,soly3,G)double(subs(f1,[t y g],[solt3 soly3 G])),solt3,soly3,G)
err31 =
1.0e-34 * Columns 1 through 10 0.0000 + 0.0000i 0.6181 + 0.4447i 0.0000 + 0.0000i 0.2703 + 0.4412i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i Column 11 -0.0000 + 0.0000i
err32 = arrayfun(@(solt3,soly3,G)double(subs(f2,[t y g],[solt3 soly3 G])),solt3,soly3,G)
err32 =
1.0e-37 * Columns 1 through 10 0.0441 + 0.0073i 0.0478 - 0.0367i 0.0698 - 0.0073i -0.0808 - 0.0661i 0.0367 + 0.0147i -0.0624 - 0.0294i 0.0220 + 0.0220i -0.1102 - 0.0220i 0.0037 + 0.0073i 0.0220 - 0.0073i Column 11 0.0184 - 0.0514i

22 Comments

Thank you @Torsten!
@Torsten Does the three splices that you created (sol1, sol2, sol3) are three solutions with t=1, t between 0.5 and 1 and t below 0.5. Which part of code takes care of this. I am new to this language and having hard time understanding it.
Torsten
Torsten on 24 Aug 2023
Edited: Torsten on 24 Aug 2023
The numerator of f2 is a polynomial of degree 3 in t and thus f2 = 0 can be solved analytically for t as a function of y and g (tsol = solve(N2==0,t,'MaxDegree',3);) . This gives three solution branches for t.
These three solution branches are then inserted in the numerator of f1 and solved for y.
It's not clear right from the beginning what these three branches will give, but it turns out that the first branch is t = 1, the second branch gives real values for t which are < 1 and the third branch gives complex values for t.
Thanks @Torsten! I just have two more questions. How to label the graph as slice 1 (for t=1), slice 2 (for t=2). And if I have to add another plot to this which has following equation -
f3 = y - 1/(1+(1/3)^(1/(1-2*g)))
(#Note:f3 =0)
Shall I just plot it
plot(G,1/(1+(1/3)^(1/(1-2*g))),'g')
Or should I vpa solve it?
Can I also label it slice 3?
Please let me know @Torsten
See above. Can you take it from here ?
You should pass the free online tutorial to learn MATLAB's basics:
Thank @Torsten. Both the solution and tutorial are super helpful
@Torsten Is there a problem with code? Why blue line (slice 2) showing bumps around 0.96? It should be continuous.
syms y t g
f1 = -2*g*t.^3/((1 - y).^(1-2*g)) + 2*g*(3 - t.^2)^2./((2 - t).*y.^(1-2*g)) + (8*g*t)./(y.^(1- 2*g)) - 2*g*(3 - 2*t).^2./((2 - t).*(1 - y).^(1-2*g));
[N1,D1] = numden(f1);
f2 = t - (((y.^(2*g))*(2/3) - (((1-y).^(2*g)).*((3 - 2*t)./(3*(2-t)))))./((y.^(2*g)).*(((3-t.^2))./(3*(2-t))) - ((1-y).^(2*g)).*(t/3)));
[N2,D2] = numden(f2);
tsol = solve(N2==0,t,'MaxDegree',3);
N1_subs_tsol1 = subs(N1,t,tsol(1));
N1_subs_tsol2 = subs(N1,t,tsol(2));
N1_subs_tsol3 = subs(N1,t,tsol(3));
G = 0.001:0.001:0.115;
soly1 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G)==0),G);
solt1 = arrayfun(@(soly1,G)subs(tsol(1),[y,g],[soly1,G]),soly1,G);
soly2 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G)==0),G);
solt2 = arrayfun(@(soly2,G)subs(tsol(2),[y,g],[soly2,G]),soly2,G);
soly3 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G)==0),G);
solt3 = arrayfun(@(soly3,G)subs(tsol(3),[y,g],[soly3,G]),soly3,G);
hold on
p1 = plot(G,soly1,'r');
%plot(G,solt1,'r')
p2 = plot(G,soly2,'b');
%plot(G,solt2,'b')
p3 = plot(G,1./(1+(1/3).^(1./(1-2*G))),'g');
legend([p1,p2,p3],[" Slice 1","Slice 2","Slice 3"])
xlabel("g")
hold off
grid on
Help @Torsten. Please
vpasolve is a numerial solver. If your equation has more than one zero, it may happen that the solver does not follow the branch of the solution it once took, but switches to another branch.
Write the commands
soly1 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G)==0),G);
solt1 = arrayfun(@(soly1,G)subs(tsol(1),[y,g],[soly1,G]),soly1,G);
soly2 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G)==0),G);
solt2 = arrayfun(@(soly2,G)subs(tsol(2),[y,g],[soly2,G]),soly2,G);
soly3 = arrayfun(@(G)vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G)==0),G);
solt3 = arrayfun(@(soly3,G)subs(tsol(3),[y,g],[soly3,G]),soly3,G);
in a loop over the elements of G and in each iteration, start "vpasolve" with an initial guess for the solution, namely the solution obtained for the previous value for G:
S = vpasolve(eqn,var,init_param)
You might be lucky and vpasolve sticks to the solution branch once taken.
syms y t g
f1 = -2*g*t.^3/((1 - y).^(1-2*g)) + 2*g*(3 - t.^2)^2./((2 - t).*y.^(1-2*g)) + (8*g*t)./(y.^(1- 2*g)) - 2*g*(3 - 2*t).^2./((2 - t).*(1 - y).^(1-2*g));
[N1,D1] = numden(f1);
f2 = t - (((y.^(2*g))*(2/3) - (((1-y).^(2*g)).*((3 - 2*t)./(3*(2-t)))))./((y.^(2*g)).*(((3-t.^2))./(3*(2-t))) - ((1-y).^(2*g)).*(t/3)));
[N2,D2] = numden(f2);
tsol = solve(N2==0,t,'MaxDegree',3);
N1_subs_tsol1 = subs(N1,t,tsol(1));
N1_subs_tsol2 = subs(N1,t,tsol(2));
N1_subs_tsol3 = subs(N1,t,tsol(3));
G = 0.01:0.001:0.12;
soly1(1) = vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G(1))==0,y);
soly2(1) = vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G(1))==0,y);
soly3(1) = vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G(1))==0,y);
for i = 2:numel(G)
soly1(i) = vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G(i))==0,y,soly1(i-1));
soly2(i) = vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G(i))==0,y,soly2(i-1));
soly3(i) = vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G(i))==0,y,soly3(i-1));
end
solt1 = arrayfun(@(soly1,G)subs(tsol(1),[y,g],[soly1,G]),soly1,G);
solt2 = arrayfun(@(soly2,G)subs(tsol(2),[y,g],[soly2,G]),soly2,G);
solt3 = arrayfun(@(soly3,G)subs(tsol(3),[y,g],[soly3,G]),soly3,G);
hold on
p1 = plot(G,soly1,'r');
%plot(G,solt1,'r')
p2 = plot(G,soly2,'b');
%plot(G,solt2,'b')
p3 = plot(G,1./(1+(1/3).^(1./(1-2*G))),'g');
legend([p1,p2,p3],[" Slice 1","Slice 2","Slice 3"])
xlabel("g")
hold off
grid on
% Check the solutions
err11 = arrayfun(@(solt1,soly1,G)double(subs(f1,[t y g],[solt1 soly1 G])),solt1,soly1,G)
err11 = 1×111
1.0e-38 * 0.0092 0.0092 0 0 0.0184 -0.0184 -0.0092 0.0092 -0.0184 0 0.0276 -0.0092 0.0092 0 -0.0092 0.0184 -0.0184 0 0.0367 0 0 0.0184 0.0367 -0.0184 0.0367 -0.0184 0 0.0184 0 -0.0551
err12 = arrayfun(@(solt1,soly1,G)double(subs(f2,[t y g],[solt1 soly1 G])),solt1,soly1,G)
err12 = 1×111
1.0e-39 * 0 0 0 0 0 0 0.1837 0 0 0 0 0 0.1837 0 0 0 0 0 0.1837 0 0.1837 0.1837 0 0 0 0.1837 0 0 0 0
err21 = arrayfun(@(solt2,soly2,G)double(subs(f1,[t y g],[solt2 soly2 G])),solt2,soly2,G)
err21 = 1×111
1.0e-37 * 0.1527 -0.0743 0.1052 -0.0202 0.1185 0.0351 -0.1309 -0.1070 0.0207 -0.0041 0.0523 0.0133 0.0225 0.0388 0.1111 0.0556 -0.0533 -0.0069 -0.0441 -0.0445 -0.0096 -0.0680 -0.0459 0.0721 -0.0046 -0.0684 -0.0354 -0.0983 0.0606 -0.0335
err22 = arrayfun(@(solt2,soly2,G)double(subs(f2,[t y g],[solt2 soly2 G])),solt2,soly2,G)
err22 = 1×111
1.0e-37 * 0.2480 -0.1102 0.1451 -0.0257 0.1378 0.0422 -0.1322 -0.1010 0.0202 0 0.0441 0.0147 0.0202 0.0276 0.0771 0.0331 -0.0349 -0.0018 -0.0202 -0.0220 -0.0055 -0.0312 -0.0220 0.0331 -0.0073 -0.0312 -0.0147 -0.0441 0.0312 -0.0147
err31 = arrayfun(@(solt3,soly3,G)double(subs(f1,[t y g],[solt3 soly3 G])),solt3,soly3,G)
err31 =
1.0e-33 * Columns 1 through 10 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i Columns 11 through 20 -0.0000 - 0.0000i -0.1191 - 0.0064i -0.0240 - 0.0034i -0.0051 - 0.0012i -0.0012 - 0.0004i -0.0003 - 0.0001i -0.0001 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i Columns 21 through 30 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i Columns 31 through 40 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i Columns 41 through 50 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i Columns 51 through 60 -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i Columns 61 through 70 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i Columns 71 through 80 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i Columns 81 through 90 -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i Columns 91 through 100 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i Columns 101 through 110 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i Column 111 0.0000 + 0.0000i
err32 = arrayfun(@(solt3,soly3,G)double(subs(f2,[t y g],[solt3 soly3 G])),solt3,soly3,G)
err32 =
1.0e-34 * Columns 1 through 10 -0.1496 + 0.1622i -0.0040 - 0.0236i 0.3863 - 0.0608i -0.0743 + 0.0123i -0.0257 - 0.0089i -0.0955 + 0.0294i 0.0026 + 0.0102i 0.0588 - 0.0052i 0.0578 + 0.0140i -0.0154 + 0.0163i Columns 11 through 20 -0.0232 + 0.0049i 0.0080 + 0.0101i 0.0125 - 0.0049i 0.0005 - 0.0018i -0.0013 - 0.0007i -0.0182 - 0.0042i -0.0147 + 0.0026i -0.0022 + 0.0030i 0.0040 - 0.0038i 0.0060 - 0.0005i Columns 21 through 30 0.0004 + 0.0028i -0.0011 - 0.0026i 0.0085 - 0.0026i 0.0026 - 0.0024i 0.0068 - 0.0004i -0.0053 + 0.0023i -0.0021 - 0.0006i -0.0008 - 0.0016i 0.0040 + 0.0015i -0.0006 - 0.0025i Columns 31 through 40 -0.0016 + 0.0005i 0.0007 - 0.0001i -0.0024 - 0.0011i -0.0010 - 0.0017i 0.0012 + 0.0003i 0.0004 - 0.0004i 0.0006 + 0.0001i -0.0008 + 0.0005i 0.0012 - 0.0003i -0.0012 + 0.0005i Columns 41 through 50 -0.0003 + 0.0006i 0.0003 + 0.0009i -0.0000 - 0.0000i -0.0000 + 0.0003i 0.0003 - 0.0003i -0.0005 - 0.0005i -0.0004 + 0.0004i -0.0001 - 0.0004i -0.0001 + 0.0003i -0.0002 - 0.0000i Columns 51 through 60 -0.0002 + 0.0001i 0.0004 + 0.0001i -0.0000 + 0.0000i 0.0002 - 0.0001i -0.0001 - 0.0002i -0.0001 + 0.0000i 0.0002 + 0.0002i 0.0001 - 0.0001i -0.0002 + 0.0000i 0.0001 + 0.0001i Columns 61 through 70 -0.0000 - 0.0001i -0.0001 - 0.0001i 0.0002 + 0.0001i -0.0002 + 0.0000i -0.0000 + 0.0002i -0.0001 - 0.0001i -0.0001 + 0.0002i -0.0000 - 0.0000i 0.0001 - 0.0002i -0.0002 - 0.0001i Columns 71 through 80 -0.0002 - 0.0000i -0.0000 - 0.0001i 0.0001 + 0.0001i 0.0001 - 0.0001i -0.0000 - 0.0000i 0.0001 + 0.0001i -0.0001 + 0.0000i 0.0001 + 0.0001i 0.0000 + 0.0001i -0.0000 + 0.0001i Columns 81 through 90 -0.0000 - 0.0000i -0.0000 - 0.0001i -0.0001 - 0.0001i -0.0000 + 0.0001i 0.0001 + 0.0000i -0.0001 + 0.0000i 0.0001 + 0.0000i -0.0000 + 0.0001i -0.0000 - 0.0000i -0.0000 - 0.0000i Columns 91 through 100 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0001 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 - 0.0000i 0.0000 - 0.0000i -0.0001 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i Columns 101 through 110 -0.0000 - 0.0001i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i Column 111 -0.0000 + 0.0000i

Thank you @Torsten. One last question. If I have to go one step ahead and want Matlab to substitute the pair of solution [that I got by solving f1 and f2 for each g in range 0.01 to 0.12) with (t=1 with associated y) and (t<1 with associated y) into following function -

f4 = (((3-t^2)^2)*(y^(2 g))+ ((3 - 2 t)^2)*(1 - k)^(2 g))/(18 (2 - t)) + (1/18) (4 t (k^(2 g)) + t^3 ((1 - k)^(2 g)))

And generate me two curves - one with solution of above function when t=1 and other when t<1 against g on x axis. How can I do it.

Help @Torsten

Torsten
Torsten on 1 Sep 2023
Edited: Torsten on 1 Sep 2023
You've got pairs (t1,y1) (t1 = 1) and (t2,y2) (t2 < 1) that simultaneously satisfy f1 and f2. Now if you insert these pairs into f4 (= 0, I suppose), you get a kind of error in how far these (t1,y1) and (t2,y2) satisfy f4. So I don't understand why you talk about "solution of above function when t=1 and other when t<1". You cannot get a "solution" from one equation in two unknowns. Or do you want to insert the t1 and t2 into f4 = 0 and solve for y ? And while answering, please correct f4 - there seem to be several (multiplication ?) signs missing.

f4 is not 0. f4 is just a function of t and y and I want the value of this function at Solution of t and y. f4 is an objective function.

   f4 = (((3-t^2)^2)*(y^(2*g))+ 
  ((3 - 2*t)^2)*(1 - 
   k)^(2*g))/(18*(2 - t)) + 
  (1/18)*(4*t*(k^(2*g)) + 
  (t^3)*((1 - k)^(2*g)))

I want to plot the value of objective function against g - two separate lines: one at t=1, other at t<1. Please help @Torsten

@Torsten Does my question make sense now? Please help
Next time I'll no longer give help if you show no effort. This question was really not too hard to manage it yourself.
syms y t g
f1 = -2*g*t.^3/((1 - y).^(1-2*g)) + 2*g*(3 - t.^2)^2./((2 - t).*y.^(1-2*g)) + (8*g*t)./(y.^(1- 2*g)) - 2*g*(3 - 2*t).^2./((2 - t).*(1 - y).^(1-2*g));
[N1,D1] = numden(f1);
f2 = t - (((y.^(2*g))*(2/3) - (((1-y).^(2*g)).*((3 - 2*t)./(3*(2-t)))))./((y.^(2*g)).*(((3-t.^2))./(3*(2-t))) - ((1-y).^(2*g)).*(t/3)));
[N2,D2] = numden(f2);
f4 = (((3-t.^2).^2).*(y.^(2*g))+ ((3 - 2*t).^2).*(1 - y).^(2*g))./(18*(2 - t)) + (1/18)*(4*t.*(y.^(2*g)) + (t.^3).*((1 - y).^(2*g)));
tsol = solve(N2==0,t,'MaxDegree',3);
N1_subs_tsol1 = subs(N1,t,tsol(1));
N1_subs_tsol2 = subs(N1,t,tsol(2));
N1_subs_tsol3 = subs(N1,t,tsol(3));
G = 0.01:0.001:0.12;
soly1(1) = vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G(1))==0,y);
soly2(1) = vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G(1))==0,y);
soly3(1) = vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G(1))==0,y);
for i = 2:numel(G)
soly1(i) = vpasolve(subs(N1_subs_tsol1/(y*(1-y)^(2*g)),g,G(i))==0,y,soly1(i-1));
soly2(i) = vpasolve(subs(N1_subs_tsol2/y/(-2*g),g,G(i))==0,y,soly2(i-1));
soly3(i) = vpasolve(subs(N1_subs_tsol3/y/(-2*g),g,G(i))==0,y,soly3(i-1));
end
solt1 = arrayfun(@(soly1,G)subs(tsol(1),[y,g],[soly1,G]),soly1,G);
solt2 = arrayfun(@(soly2,G)subs(tsol(2),[y,g],[soly2,G]),soly2,G);
solt3 = arrayfun(@(soly3,G)subs(tsol(3),[y,g],[soly3,G]),soly3,G);
f41 = arrayfun(@(solt1,soly1,G)subs(f4,[t,y,g],[solt1,soly1,G]),solt1,soly1,G);
f42 = arrayfun(@(solt2,soly2,G)subs(f4,[t,y,g],[solt2,soly2,G]),solt2,soly2,G);
hold on
p1 = plot(G,soly1,'r');
%plot(G,solt1,'r')
p2 = plot(G,soly2,'b');
%plot(G,solt2,'b')
p3 = plot(G,1./(1+(1/3).^(1./(1-2*G))),'g');
p41 = plot(G,f41,'--r');
p42 = plot(G,f42,'--b');
legend([p1,p41,p2,p42,p3],["Slice 1","Slice 1' ","Slice 2","Slice 2' ","Slice 3"],'Position',[0.7 0.4 0.1 0.2])
xlabel("g")
hold off
grid on
% Check the solutions
err11 = arrayfun(@(solt1,soly1,G)double(subs(f1,[t y g],[solt1 soly1 G])),solt1,soly1,G)
err11 = 1×111
1.0e-38 * 0.0092 0.0092 0 0 0.0184 -0.0184 -0.0092 0.0092 -0.0184 0 0.0276 -0.0092 0.0092 0 -0.0092 0.0184 -0.0184 0 0.0367 0 0 0.0184 0.0367 -0.0184 0.0367 -0.0184 0 0.0184 0 -0.0551
err12 = arrayfun(@(solt1,soly1,G)double(subs(f2,[t y g],[solt1 soly1 G])),solt1,soly1,G)
err12 = 1×111
1.0e-39 * 0 0 0 0 0 0 0.1837 0 0 0 0 0 0.1837 0 0 0 0 0 0.1837 0 0.1837 0.1837 0 0 0 0.1837 0 0 0 0
err21 = arrayfun(@(solt2,soly2,G)double(subs(f1,[t y g],[solt2 soly2 G])),solt2,soly2,G)
err21 = 1×111
1.0e-37 * 0.1527 -0.0743 0.1052 -0.0202 0.1185 0.0351 -0.1309 -0.1070 0.0207 -0.0041 0.0523 0.0133 0.0225 0.0388 0.1111 0.0556 -0.0533 -0.0069 -0.0441 -0.0445 -0.0096 -0.0680 -0.0459 0.0721 -0.0046 -0.0684 -0.0354 -0.0983 0.0606 -0.0335
err22 = arrayfun(@(solt2,soly2,G)double(subs(f2,[t y g],[solt2 soly2 G])),solt2,soly2,G)
err22 = 1×111
1.0e-37 * 0.2480 -0.1102 0.1451 -0.0257 0.1378 0.0422 -0.1322 -0.1010 0.0202 0 0.0441 0.0147 0.0202 0.0276 0.0771 0.0331 -0.0349 -0.0018 -0.0202 -0.0220 -0.0055 -0.0312 -0.0220 0.0331 -0.0073 -0.0312 -0.0147 -0.0441 0.0312 -0.0147
err31 = arrayfun(@(solt3,soly3,G)double(subs(f1,[t y g],[solt3 soly3 G])),solt3,soly3,G)
err31 =
1.0e-33 * Columns 1 through 10 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i Columns 11 through 20 -0.0000 - 0.0000i -0.1191 - 0.0064i -0.0240 - 0.0034i -0.0051 - 0.0012i -0.0012 - 0.0004i -0.0003 - 0.0001i -0.0001 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i Columns 21 through 30 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i Columns 31 through 40 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i Columns 41 through 50 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i Columns 51 through 60 -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i Columns 61 through 70 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i Columns 71 through 80 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i Columns 81 through 90 -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i Columns 91 through 100 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i Columns 101 through 110 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i Column 111 0.0000 + 0.0000i
err32 = arrayfun(@(solt3,soly3,G)double(subs(f2,[t y g],[solt3 soly3 G])),solt3,soly3,G)
err32 =
1.0e-34 * Columns 1 through 10 -0.1496 + 0.1622i -0.0040 - 0.0236i 0.3863 - 0.0608i -0.0743 + 0.0123i -0.0257 - 0.0089i -0.0955 + 0.0294i 0.0026 + 0.0102i 0.0588 - 0.0052i 0.0578 + 0.0140i -0.0154 + 0.0163i Columns 11 through 20 -0.0232 + 0.0049i 0.0080 + 0.0101i 0.0125 - 0.0049i 0.0005 - 0.0018i -0.0013 - 0.0007i -0.0182 - 0.0042i -0.0147 + 0.0026i -0.0022 + 0.0030i 0.0040 - 0.0038i 0.0060 - 0.0005i Columns 21 through 30 0.0004 + 0.0028i -0.0011 - 0.0026i 0.0085 - 0.0026i 0.0026 - 0.0024i 0.0068 - 0.0004i -0.0053 + 0.0023i -0.0021 - 0.0006i -0.0008 - 0.0016i 0.0040 + 0.0015i -0.0006 - 0.0025i Columns 31 through 40 -0.0016 + 0.0005i 0.0007 - 0.0001i -0.0024 - 0.0011i -0.0010 - 0.0017i 0.0012 + 0.0003i 0.0004 - 0.0004i 0.0006 + 0.0001i -0.0008 + 0.0005i 0.0012 - 0.0003i -0.0012 + 0.0005i Columns 41 through 50 -0.0003 + 0.0006i 0.0003 + 0.0009i -0.0000 - 0.0000i -0.0000 + 0.0003i 0.0003 - 0.0003i -0.0005 - 0.0005i -0.0004 + 0.0004i -0.0001 - 0.0004i -0.0001 + 0.0003i -0.0002 - 0.0000i Columns 51 through 60 -0.0002 + 0.0001i 0.0004 + 0.0001i -0.0000 + 0.0000i 0.0002 - 0.0001i -0.0001 - 0.0002i -0.0001 + 0.0000i 0.0002 + 0.0002i 0.0001 - 0.0001i -0.0002 + 0.0000i 0.0001 + 0.0001i Columns 61 through 70 -0.0000 - 0.0001i -0.0001 - 0.0001i 0.0002 + 0.0001i -0.0002 + 0.0000i -0.0000 + 0.0002i -0.0001 - 0.0001i -0.0001 + 0.0002i -0.0000 - 0.0000i 0.0001 - 0.0002i -0.0002 - 0.0001i Columns 71 through 80 -0.0002 - 0.0000i -0.0000 - 0.0001i 0.0001 + 0.0001i 0.0001 - 0.0001i -0.0000 - 0.0000i 0.0001 + 0.0001i -0.0001 + 0.0000i 0.0001 + 0.0001i 0.0000 + 0.0001i -0.0000 + 0.0001i Columns 81 through 90 -0.0000 - 0.0000i -0.0000 - 0.0001i -0.0001 - 0.0001i -0.0000 + 0.0001i 0.0001 + 0.0000i -0.0001 + 0.0000i 0.0001 + 0.0000i -0.0000 + 0.0001i -0.0000 - 0.0000i -0.0000 - 0.0000i Columns 91 through 100 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0001 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 - 0.0000i 0.0000 - 0.0000i -0.0001 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i Columns 101 through 110 -0.0000 - 0.0001i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0001 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i Column 111 -0.0000 + 0.0000i

Why I am getting this error -

   Error using arrayfun
All of the input arguments must be of the same
size and shape.
Previous inputs had size 118 in dimension 2.
Input #3 has size 111

And I am really sorry - there is no k in f4, it is y. It is my typo.

   f4 = (((3-t.^2).^2).*(y.^(2*g))+ ((3 - 2*t).^2).*(1 - y).^(2*g))./(18*(2 - t)) + (1/18)*(4*t.*(y.^(2*g)) + (t.^3).*((1 - y).^(2*g)));

Please help @Torsten

Torsten
Torsten on 1 Sep 2023
Edited: Torsten on 1 Sep 2023
I don't have such a problem (see above).
Maybe you chose a different G and didn't get a solution y for all of your G-values.
@Torsten I just copied your code as it is and I am still getting the error. Is it because I am running on mobile?
Can you please run it for me? Give me graph of just p41 and p42 and value of g from 0.01 to 0.49? Or at least g from 0.01 to 0.15
Sorry for disturbing you. Thank you @Torsten for everything.
Help @Torsten Please
It works up to g = 0.12. For higher values of g, vpasolve has trouble to find a solution. I don't have the time to analyze for the reasons, sorry.
@Torsten Can you please plot me till 0.12? Plot of just slice 1' abs slice 2' because I still can't run it.
You can edit the code according to your needs using the edit pencil and run it here with MATLAB online using the green RUN arrow.
@Torsten Thank you a ton

Sign in to comment.

More Answers (1)

Vikas
Vikas on 23 Aug 2023
I read you question on the above that I understand you want to plot one to many graph in one graph.
For this you use "hold on" command
This is use to hold the previous graph after plot and then plot another graph that affect in same graph.
example:-
plot(3,6)
hold on
plot(44,23)
This two plot draw a graph in same only one graph.
As I wish you understand this answer?.

Products

Release

R2023a

Community Treasure Hunt

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

Start Hunting!