Legend for varying parameters in a loop

Is it possible to have legends such that for xi=-0.06 (I have 4 different plots for various times). I want the lines to have the same color(say blue) and for xi=-0.05 (I have 4 different plots for various times) I want the lines to have the same color(say red)? See the lines concern and I have also attached my codes

f1 = @(xi) theta(1:round(nt/4):nt,:);

hold on

for i = 1:numel(xi)

plot(z, f1(xi(i)))

xlabel('t (seconds)')

ylabel('\theta(z,t)(rad)')

legend('\xi=-0.06', '\xi=-0.05')

end

hold off

Answer 1

[Walter Roberson]

plot(z, f1(xi(i)))

assign the result of f1(xi(i)) to a variable, say f1i. Then

h(i, :)  = plot(z, f1i(:, 1), 'b', f1i(:, 2), 'r');

Now delete the legend call inside the loop.

After the loop

legend(h(1,:), {'\xi=-0.06', '\xi=-0.05'})

[University Glasgow]

Thank you.

f1 = @(xi) theta(1:round(nt/4):nt,:);

hold on

for i = 1:numel(xi)

f1i = f1(xi(i));

h(i, :) = plot(z, f1i(:, 1), 'b', z, f1i(:, 2), 'r');

xlabel('t (seconds)')

ylabel('\theta(z,t)(rad)')

end

legend(h(1,:), {'\xi=-0.06', '\xi=-0.05'})

hold off

I tried it but I received this error:

Error using plot Vectors must be the same length.

Error in Revised_Case1_Ijuptilk (line 93)

h(i, :) = plot(z, f1i(:, 1), 'b', f1i(:, 2), 'r');

[Walter Roberson]

your code is different than your error message. the code in your error message does not have the second z

I just noticed that your f1 ignores its input so you are drawing the same thing each iteration, which sounds wrong.

how many columns does theta have? How long is z?

[Torsten]

f1 = @(xi) theta(1:round(nt/4):nt,:);

Your f1 does not depende on xi. So you'll always get back the same part of the theta matrix.

Concerning the error message: before the plot command, insert the lines

size(z)
size(f1i)

What do you get back from MATLAB ?

[University Glasgow]

The size(z) is 31 and f1i is 5 by 31

[University Glasgow]

Still same thing without the second z.

f1 = @(xi) theta(1:round(nt/4):nt,:);

hold on

for i = 1:numel(xi)

f1i = f1(xi(i));

h(i, :) = plot(z, f1i(:, 1), 'b', f1i(:, 2), 'r');

xlabel('t (seconds)')

ylabel('\theta(z,t)(rad)')

end

legend(h(1,:), {'\xi=-0.06', '\xi=-0.05'})

hold off

[Torsten]

size(z) cannot be 31. It must either be 31x1 or 1x31, and which of the two is correct is important.

[University Glasgow]

Sorry, size(z) is 1x31

[Torsten]

h(i, :)  = plot(z, f1i(1,:), 'b', f1i(2,:), 'r');

instead of

h(i, :)  = plot(z, f1i(:, 1), 'b', f1i(:, 2), 'r');

[University Glasgow]

Thank you I just did now but is working for only plot(z, f1i(:, 1))

[Walter Roberson]

h(i, :)  = plot(z, f1i(1,:), 'b', z, f1i(2,:), 'r');

[Torsten]

... and as the f1i are all the same, colors won't matter :-)

[University Glasgow]

With the code below, I don't think f1i are the same:

f1 = @(xi) theta(1:round(nt/4):nt,:);

hold on

for i = 1:numel(xi)

plot(z, f1(xi(i)))

xlabel('t (seconds)')

ylabel('\theta(z,t)(rad)')

legend('\xi=-0.06', '\xi=-0.05')

end

hold off

You can see from the figures I attached: The figure for xi =[-0.06, -0.05] and xi =[-0.06, -0.06] are different. Remember I have loop over my solver for xi.

[Torsten]

How should the f1i be different if the function f1 does not depend on xi ?

This can only happen if theta has changed because of an outer loop over the xi.

But then, your inner loop from above - also over xi - cannot make sense.

Please include the code you use at the moment.

[Torsten]

%% initialization
clear 
clc
close all


%% Model parameters
global N Phi A B C G h
k1 = 6*10^(-12);            %  elastic constant
eta1 = 0.0240;              %  viscosity
alpha3 = -0.001104;         %  viscosity
gamma1 = 0.1093;            %  viscosity
Theta = 0.0001;
d = 0.0002;

% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
N = 30; % N must be even


%% Numerical setup

% step size
h = d/N;
% t span
tspan = [0:100];
nt=length(tspan);

% range of z
z=linspace(0,d,N+1);

% initial conditions
theta0 = Theta*sin(pi*z/d);
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];

% Constant mass matrix M: We use M to represent the left hand side of the system of equations.
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];

% Boundary Conditions
Phi = 0;


%% ode solver
for xi = [-0.06, -0.05]
    options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
    [t,y] = ode15s(@(t,y)lcode1(t, y, xi),tspan,u0, options);
    
    % Extract the solution for theta and v
    theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
    v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
    
    
    
    % Extract theta solution for different value
    mytheta= xlsread('thetasolxi1.xlsx');
    thetasolxi = mytheta(:,[2, 4, 6]);
    z1 = mytheta(:,1);
    
    % Extract v solution for different value
    myvol= xlsread('vsolxi1.xlsx');
    vsolxi = myvol(:,[2, 4, 6]);
    z2 = myvol(:,1);
    
    % theta at the middle of the layer (i.e., z=d/2)
    theta_middle = theta(:, N/2);     
    
    
    %% Extract theta_middel data from maple
    mydata= xlsread('DatafromMaple_middle.xlsx');
    theta_middlemaple = mydata(:,2);
    time = mydata(:,1);
    
    
    % Extract r(q) data from maple
    rqdata= xlsread('rqdata.xlsx');
    rq= rqdata(:,2);
    qsize = rqdata(:,1);
    
    % % theta plot for different xi values
    f1 = @(xi) theta(1:round(nt/4):nt,:);
    hold on
    for i = 1:numel(xi)
        plot(z, f1(xi(i)))
        xlabel('t (seconds)')
        ylabel('\theta(z,t)(rad)')
        legend('\xi=-0.06', '\xi=-0.05')
    end
    hold off
end

% %% Plot the solution. 
% figure
% subplot(2,1,1)
% plot(z, theta(1:round(nt/10):nt,:))
% xlabel('z(\mum)')
% ylabel('\theta(z,t)(rad)')
% 
% subplot(2,1,2)
% plot(z, v(1:round(nt/10):nt,:))
% xlabel('z(\mum)')
% ylabel('v(z,t)(m/s)')
% %legend('t_0 =0', 't_1=10', 't_2=20', 't_3=30', 't_4=40', 'Location','bestoutside')
% 
% % matlab and maple for thetamiddle plot
% figure
% hold on
% plot(t,theta_middle, 'r', MarkerSize=2)
% plot(time,theta_middlemaple, 'b--', MarkerSize=2)
% xlabel('t (seconds)')
% ylabel('\theta(d/2,t)(rad)')
% 
% 
% % matlab and maple for theta plot for different xi values
% figure
% hold on
% plot(z, theta(end, :))
% plot(z1, thetasolxi./2, 'r', MarkerSize=2)
% xlabel('t (seconds)')
% ylabel('\theta(z,t)(m/s)')
% 
% 
% % matlab and maple for v plot for different xi values
% figure
% hold on
% plot(z, v(end, :))
% plot(z2, vsolxi./2, 'r', MarkerSize=2)
% xlabel('t (seconds)')
% ylabel('v(z,t)(m/s)')
% 
% 
% % r(q) equation
% alpha=1-alpha3^2/gamma1/eta1;
% q=0:0.001:(5*pi);
% r=q-(1-alpha)*tan(q)+(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1).*tan(q);
% 
% % Plotting r(q)
% figure
% plot(q,r, 'b--', MarkerSize=2)
% xlabel('t (seconds)')
% ylabel('\theta(d/2,t)(rad)')
% %legend('\xi=-0.02','\xi=0.01', '\xi=0.02')
% axis([0 5*pi -20 20])

%% Definition of Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t, y ,xi)
global N Phi A B C G h

% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));

% theta equations
% for the positon to the right of the left hand boundary z=0
rhsode(1,1) =     (A/(h^2))*(theta(2)-2*theta(1)+ Phi)          - (B/(2*h))*(v(2)-0);
% for all other internal positions 0<z<d
for i=2:(N-2)
    rhsode(i,1) = (A/(h^2))*(theta(i+1)-2*theta(i)+theta(i-1))  - (B/(2*h))*(v(i+1)-v(i-1)); 
end 
% for the positon to the left of the right hand boundary z=d
rhsode(N-1,1)   = (A/(h^2))*(Phi-2*theta(N-1)+ theta(N-2))      - (B/(2*h))*(0-v(N-2));

% v equations [REMEMBER theta the RHS index for the v equations is N-1 more than the indices of the theta and v variables in this function]
% for the two positons to the right of the left hand boundary z=0
rhsode(N,1)     = (G/(h^3))*(-Phi + 3*theta(1) - 3*theta(2) + theta(3))            + (C/(h^2))*(v(2) -2*v(1)+ 0)             + (xi/(2*h))*(theta(2)-Phi);
rhsode(N+1,1)   = (G/(2*h^3))*(Phi - 2*theta(1) + 2*theta(3)- theta(4))            + (C/(h^2))*(v(3) -2*v(2) + v(1))         + (xi/(2*h))*(theta(3)-theta(1));
% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
    rhsode(i,1) = (G/(2*h^3))*(theta(i-2-(N-1)) - 2*theta(i-1-(N-1)) + 2*theta(i+1-(N-1))- theta(i+2-(N-1)))   +(C/(h^2))*(v(i+1-(N-1)) -2*v(i-(N-1)) + v(i-1-(N-1)))          + (xi/(2*h))*(theta(i+1-(N-1))-theta(i-1-(N-1)));
end
% for the two positons to the left of the right hand boundary z=d
rhsode(2*N-3,1) = (G/(2*h^3))*(theta(N-4) - 2*theta(N-3) + 2*theta(N-1) - Phi)     +(C/(h^2))*(v(N-1) -2*v(N-2) + v(N-3))    + (xi/(2*h))*(theta(N-1)-theta(N-3));
rhsode(2*N-2,1) = (G/(h^3))*(-theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi)      +(C/(h^2))*(0 -2*v(N-1) + v(N-2))         + (xi/(2*h))*(Phi-theta(N-2));      

end

[Torsten]

It makes no sense to make a plot loop over the xi before the calculations for all the xi have finished .

Save the solution for theta and v first in 3d-arrays (the 3rd dimension respects the different xi) and plot after all calculations have finished outside the calculation loop

for xi = [-0.06, -0.05]
    ...
end

[University Glasgow]

Please how I theta and v first in 3d-arrays?

[Torsten]

counter = 0;
for xi = [-0.06, -0.05]
    options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
    [t,y] = ode15s(@(t,y)lcode1(t, y, xi),tspan,u0, options);
    
    % Extract the solution for theta and v
    theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
    v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
    
    % Extract theta solution for different value
    mytheta= xlsread('thetasolxi1.xlsx');
    thetasolxi = mytheta(:,[2, 4, 6]);
    z1 = mytheta(:,1);
    
    % Extract v solution for different value
    myvol= xlsread('vsolxi1.xlsx');
    vsolxi = myvol(:,[2, 4, 6]);
    z2 = myvol(:,1);
    
    % theta at the middle of the layer (i.e., z=d/2)
    theta_middle = theta(:, N/2);     
    
    
    %% Extract theta_middel data from maple
    mydata= xlsread('DatafromMaple_middle.xlsx');
    theta_middlemaple = mydata(:,2);
    time = mydata(:,1);
    
    
    % Extract r(q) data from maple
    rqdata= xlsread('rqdata.xlsx');
    rq= rqdata(:,2);
    qsize = rqdata(:,1);
    
    counter = counter + 1;
    Theta(counter,:,:) = theta;
    V(counter,:,:) = v;
end

[University Glasgow]

Wow, I see what you mean now. I really appreciate your help. Now the dimension of Theta and V is now 2x101x31 respectively. Please how do I extract Theta and V in 2D ?

[Torsten]

These are all basic questions, and you already work with MATLAB now for several weeks.

If you want to continue working with this tool, I suggest the MATLAB Onramp course:

https://de.mathworks.com/learn/tutorials/matlab-onramp.html

for ixi = 1:2
    theta = Theta(ixi,:,:)
    v = V(ixi,:,:)
    ...
end

[University Glasgow]

Thank you so much. I will definietly enrol in the course?

[University Glasgow]

Hi, pleae why does this gives me the same plot?

  xi = [-0.06, -0.1];  
    figure
    hold on
    for ixi = 1:numel(xi)
        plot(z, squeeze(THETA(ixi,1:round(nt/4):nt,:))); 
    end
    hold off

I was expecting my first plot to corresponds to xi =-0.06 and the second plot to corresponds to xi = -0.1 but I got the same plot.

[Torsten]

I don't know. I get all curves stored in the plot.

f1 = @(x) x.^2;
f2 = @(x) exp(x);
z = 0:0.01:1;
THETA(1,:,:) = [f1(z);f2(z)];
THETA(2,:,:) = 3*THETA(1,:,:);  
size(THETA)

ans = 1×3

2 2 101

figure
hold on
for ixi = 1:2
    if ixi == 1
      plot(z, squeeze(THETA(ixi,:,:)),'r'); 
    else
      plot(z, squeeze(THETA(ixi,:,:)),'b');
    end
end
grid on
hold off

[University Glasgow]

Thank you so much

[University Glasgow]

The code I have here gives legend with same color for corresponding values of xi, which is okay. But I want the legend to appear such that for xi=-0.06 (r), for xi=-0.2 (b) and for xi=-0.06 (c) with only one line representing the 3 lines instead of 3 lines for each value of xi. I don't know if this is possible.

hold on

for ixi = 1:3

if ixi == 1

plot(z, squeeze(THETA(ixi,1:round(nt/3):nt,:)),'r--', 'DisplayName','\xi=-0.06');

elseif ixi==2

plot(z, squeeze(THETA(ixi,1:round(nt/3):nt,:)),'b', 'DisplayName','\xi=-0.2');

else

plot(z, squeeze(THETA(ixi,1:round(nt/3):nt,:)),'c', 'DisplayName','\xi=-0.1');

end

xlabel('t (seconds)')

ylabel('\theta(z,t)(rad)')

end