# plot function of two variables as one variable

3 views (last 30 days)
University Glasgow on 17 Aug 2022
Commented: Torsten on 17 Aug 2022
I want to plot theta in the middle of the layer by setting z=d/2, for t=0..100. I have attached the snapshot of what I'm expecting. See my code below: %% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.06; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
d = 0.0002;
N = 20;
%% Numerical setup
% step size
h = d/N;
% t span
tspan = [0 100];
% range of z
%z=linspace(0,d/2,N+1)
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
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,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)];
%% Plot the solution.
figure
plot(z, theta(end,:))
xlabel('z(\mum)')
title('Director angle plot')
%% Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));
% theta equations
% 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 on 17 Aug 2022
Edited: Torsten on 17 Aug 2022
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.06; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
d = 0.0002;
N = 20;
%% Numerical setup
% step size
h = d/N;
% t span
tspan = [0 230];
% range of z
%z=linspace(0,d/2,N+1)
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
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,u0, options);
theta_middle = y(:,N/2);
plot(t,theta_middle) % 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)];
%% Plot the solution.
%figure
%plot(z, theta(end,:))
%xlabel('z(\mum)')
%title('Director angle plot')
%% Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));
% theta equations
% 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 on 17 Aug 2022
Sorry, but I have no experience in optimizing graphics output.
Maybe you should open a new question.

R2022a

### Community Treasure Hunt

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

Start Hunting!