plot function of two variables as one variable
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
University Glasgow
le 17 Août 2022
Commenté : Torsten
le 17 Août 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)')
ylabel('\theta(z,t)(rad)')
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
0 commentaires
Réponse acceptée
Torsten
le 17 Août 2022
Modifié(e) : Torsten
le 17 Août 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)')
%ylabel('\theta(z,t)(rad)')
%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
5 commentaires
Torsten
le 17 Août 2022
Sorry, but I have no experience in optimizing graphics output.
Maybe you should open a new question.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur General Physics dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!