Problem with for loop plotting a graph

Hi I'm really new to Matlab so you'll have to bare with me. I'm trying to plot a graph using the for function but to no avail.

I'm using the following code:

>> for a=2.0:-0.02:0.0

p=[f(a),g(a),h(a),r(a),t(a)] (where f,g,h,r,t are complicated functions of a I don't fancy writing down.)

R=roots(p)

C=imag(R)

c=C(1,1)

....

end

I'm looking to plot the complex part of the first root of each polynomial against each value of a or if its possible plotting the complex part of every root against a for every a on the same graph with different colours representing each root.

Thanks, sorry if the problem is ill described.

[madhan ravi]

attach the script file and the relevent function files

[Alexander Kimbley]

Thanks, this is the entire polynomial.

Answer 1

[Star Strider]

The Symbolic Math Toolbox is inefficient for this sort of problem.

Try this:

U = 1;
B = 0.1;
H = 0.5;
pfcn = @(a) [2*exp(2*a*(1-H)).*(1 - tanh(a*H)), -6*U*(1 + tanh(a*H)).*exp(2*a*(1 - H)) + 2*U*(1 - tanh(a*H)),(7*U^2 *(1 + tanh(a*H)) - 2*B^2 *tanh(a*H)).*exp(2*a*(1 - H)) - 5*U^2 *(1 - tanh(a*H)), exp(2*a*(1 - H)).*(2*U*B^2 *tanh(a*H) - 4*U^3 *(1 + tanh(a*H))) - (2*U*B^2 *tanh(a*H) - 4*U^3 *(1 + tanh(a*H))), U^2 *(U^2 *(1 + tanh(a*H)) - B^2 *tanh(a*H)).*exp(2*a*(1 - H)) - U^2 *(U^2 *(1 - tanh(a*H)) - B^2 *tanh(a*H))];

a=2.0:-0.02:0.0;
for k1 = 1:numel(a)
    p = pfcn(a);
    R = roots(p);
    C(:,k1) = imag(R);
end

figure
plot(a, C, '.')
grid

figure
meshc(C)
grid on

That should tell you all you need to know about the imaginary roots of your function at each ‘a’ value.

[Alexander Kimbley]

Thanks a lot! Unfortunatley I get the following error message:

"Unable to preform assignment because the size of the left side is 4 by 1 and the size of the right side is 504-1"

The problem seems to be in the "for" code.

[Star Strider]

My pleasure!

I made one small change to the ‘p’ assignment, and forgot to update my posted code:

for k1 = 1:numel(a)
    p = pfcn(a(k1));
    R = roots(p);
    C(:,k1) = imag(R);
end

I thought I’d copied the correct version. My apologies for not checking before I clicked on ‘Submit’.

[Alexander Kimbley]

Oh okay I see! Don't worry about it, it's working perffectly now! Thank you very much, it's really helped a lot!

[Star Strider]

As always, my pleasure!

[Alexander Kimbley]

Hi, this is just a follow up question but would it be possible to use a polynomail that is in factorised form. It's 6th order and obviously would be absolutley horrible to get it in the form ax^6 .....=0.

Thanks,

Alexander.

[Star Strider]

As always, my pleasure.

It depends on what you want to do with it.

A 6th-degree polynomial would not have an analytic expression for the roots, but you could certainly calculate them numerically. If it’s in a form that could be put into an anonymous function such as ‘pfcn’, it shouldn’t be much of a problem. By the same token, you can always use the conv (link) function to construct it from its factors.

[Alexander Kimbley]

Thanks for the quick reply! Well I need to use the same code you provided for me however pfcn is now of 6th order. I'll attach it, let me know what you think.

[Star Strider]

Is it a function of ‘a’, ‘x’, or one of the others? What would typical values be for the parameters, and the range of the independent variable?

[Alexander Kimbley]

X Is a complex number, B constant, a>0 constant but I vary in the analysis after (hence the loop) and H is also constant. The actual value of the function isnt needed but the coeficients of x^n are to be used in our pfcn function.

The reason why I need H and B to remain is so I can change them to see how the magnitude of the complex part of x changes with them.

[Star Strider]

There should be six roots, but there is only one. You appear to be doing a straghtforward multiplication rather than doing a colvolution:

B = rand;

a = rand;

H = rand;

pfcn = @(x) [B^2+exp(2*a*H)*(B^2-2*x.^2)].*[B^2-(1-x).^2].*[((1-x).^2)*(1+tanh(a*(1-H)))-B^2-a*(((1-x).^2).*(1-tanh(a*(1-H)))-B^2)]-(B^2-2*x.^2)*(B^2-a*(B^2-2*x.^2)).*[((1-x).^2)*(1+tanh(a*(1-H)))-B^2+exp(-2*a*H)*(((1-x).^2)*(1-tanh(a*(1-H)))-B^2)];

rts = fsolve(pfcn, (rand+1i*rand));

You need to convolve the various factors to produce the polynomial you want. You likely only need to do this once, and the Symbolic Math Toolbox may be the best approach. One option would be:

syms x a b c 
r1 = x + a;
r2 = x + b;
r3 = x + c;
p(x) = expand(r1 * r2 * r3);
p = collect(p, x)

producing:

p(x) =
x^3 + (a + b + c)*x^2 + (a*b + a*c + b*c)*x + a*b*c

Using it with the matlabFunction function would produce an anonymous function for ‘pfcn’:

pfcn = matlabFunction(p)

producing:

pfcn =

  function_handle with value:

    @(x,a,b,c)x.^3+(a+b+c).*x.^2+(a.*b+a.*c+b.*c).*x+a.*b.*c

or more simply (after a bit of editing):

    pfcn = @(x,a,b,c)x.^3+(a+b+c).*x.^2+(a.*b+a.*c+b.*c).*x+a.*b.*c;

The coefficients would not need to be expressly used as part of the argument list if you wanted ‘pfcn’ to inherit them from your workspace.

[Alexander Kimbley]

Right, thanks. I've manged to ge tthe following expresson for pfcn and subbed into our code however there is an error. Let me know what you think. I'll attach the code to this.

Thanks.

[Star Strider]

As always, my pleasure.

This now works, and produces interesting plots:

B=0.5;
H=0.2;
pfcn = @(a) [((4*exp(2*H*a) - 2*exp(4*H*a) - (2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (2*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (4*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (4*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + 2) + (8*exp(4*H*a) - 12*exp(2*H*a) + (4*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (4*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (12*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (8*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (12*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (8*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - 4));
    + ((14*exp(2*H*a) - 12*exp(4*H*a) - (2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + B^2*exp(-2*H*a) - 9*B^2*exp(2*H*a) + 5*B^2*exp(4*H*a) - 5*B^2 + (2*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (14*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (12*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (3*B^2*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (14*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (12*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (3*B^2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (B^2*exp(-2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (7*B^2*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (3*B^2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (B^2*exp(-2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (7*B^2*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (3*B^2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + 2));
    + ((8*exp(4*H*a) - 8*exp(2*H*a) - 2*B^2*exp(-2*H*a) + 14*B^2*exp(2*H*a) - 12*B^2*exp(4*H*a) + 8*B^2 + (8*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (8*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (4*B^2*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (8*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (8*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (4*B^2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (2*B^2*exp(-2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (18*B^2*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (8*B^2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (2*B^2*exp(-2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (18*B^2*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (8*B^2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1)));
    + ((2*exp(2*H*a) - 2*exp(4*H*a) + B^2*exp(-2*H*a) - 7*B^2*exp(2*H*a) + 10*B^2*exp(4*H*a) - 2*B^4*exp(-2*H*a) + 6*B^4*exp(2*H*a) - 4*B^4*exp(4*H*a) - 8*B^2 + 4*B^4 - (2*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (2*B^2*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (B^4*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (2*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (2*B^2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (B^4*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (B^2*exp(-2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (17*B^2*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (8*B^2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (B^4*exp(-2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (3*B^4*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (B^4*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (B^2*exp(-2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (17*B^2*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (8*B^2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (B^4*exp(-2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (3*B^4*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (B^4*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1)));
    + ((2*B^4*exp(-2*H*a) - 4*B^2*exp(4*H*a) - 2*B^4*exp(2*H*a) + 4*B^4*exp(4*H*a) + 4*B^2 - 4*B^4 - (4*B^2*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (2*B^4*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (4*B^2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (2*B^4*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (8*B^2*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (4*B^2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (2*B^4*exp(-2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (6*B^4*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (2*B^4*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (8*B^2*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (4*B^2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (2*B^4*exp(-2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (6*B^4*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (2*B^4*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1)));
    + (B^2*exp(4*H*a) - B^4*exp(-2*H*a) + B^4*exp(2*H*a) - 2*B^4*exp(4*H*a) + B^6*exp(-2*H*a) - B^6*exp(2*H*a) + B^6*exp(4*H*a) - B^2 + 2*B^4 - B^6 + (B^2*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (B^4*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (B^2*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (B^4*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (2*B^2*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) - (B^2*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (B^4*exp(-2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (3*B^4*exp(2*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (B^4*exp(4*H*a)*tanh(H*a))/(tanh(H*a)*tanh(a) - 1) + (2*B^2*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) + (B^2*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (B^4*exp(-2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (3*B^4*exp(2*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1) - (B^4*exp(4*H*a)*tanh(a))/(tanh(H*a)*tanh(a) - 1))];


a=10.0:-0.001:0.0;
for k1 = 1:numel(a)
    p = pfcn(a(k1));
    R = roots(p);
    C(:,k1) = imag(R);
end
figure
plot(a, C, '.')
grid
figure
meshc(C)
grid on

The problem was that you had some commas and periods in what appear to be inappropriate places, and also concatenation errors created by the spaces between the terms and operators in each line. (The solution for the concatenation errors is to remove the commas and periods, and for the spaces, to put parentheses around each complete line.)

[Alexander Kimbley]

Thanks a lot! You've been such a godsend!

[Star Strider]

As always, my pleasure! Thank you!