Extract of data for the last time step

3 次查看（过去 30 天）
University Glasgow
University Glasgow 2022-12-2
评论： University Glasgow 2022-12-2
I have tried to extract theta at the middle of the array for all t. Please, how can I extract theta middle that correspond to the last time step? See the my code:
N = 100;
%% Numerical setup
% step size
h = d/N;
tspan = 0:2000; % Nonlinear case
nt=length(tspan);
% range of z
z=linspace(0,d,N+1);
% initial conditions
Theta = 0.0001;
theta0 = Theta*sin(pi*z/d); % Initial condition: 0
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];
% Matrix M
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];
% Boundary Conditions
Phi_b = 0;
%% ode solver
counter = 0;
Xis = xi;
for xiidx = 1:length(Xis)
xi = Xis(xiidx);
for Uis = u
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
tic
[t,y] = ode15s(@(t,y)lcode1(t, y), tspan,u0, options);
toc
% Extract the solution for theta and v
theta = [Phi_b*ones(length(t), 1) y(:,1:N-1) Phi_b*ones(length(t), 1)];
v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
% Extract theta and v data for various xi
counter = counter + 1;
THETA(counter,:,:) = theta;
V(counter,:,:) = v;
% theta at the middle of the layer (i.e., z=d/2)
theta_middle(counter,:,:) = theta(:, N/2 +1);
end
end
  2 个评论
Jan
Jan 2022-12-2
I struggle with the term "theta middle at the last for last time step". What does it mean?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

William Rose
William Rose 2022-12-2
@University Glasgow, Please post the simplest possible code example that illustrates the issue you have. When you include code, you can include it as actual code, by selecting it and then clicking the "Code" button at the top of the editing window. When you include a function, please also include a (short) script illustrating a call to that function.
You have two nested for loops. Inside the inner loop, you call ode45(). ode45() returns column vector t and array y. Array theta is the same as array y, expcept a column has been added on the right of y and on the left of y. (The added columns have the constant value Phi_b.) You create a 3D array THETA. Each slice of THETA contains the 2D array theta from one pass through the inner loop. THETA is never used, so I'm not sure why it is there. Since tspan is a vector of length 2001 on every pass, t and y will have length 2001 on each pass.
Your code includes the lines
theta_middle(counter,:,:) = theta(:, N/2 +1);
theta_middle12(counter,:,:) = theta(:, N/2 +1);
theta_middle13(counter,:,:) = theta(:, N/2 +1);
theta_middle14(counter,:,:) = theta(:, N/2 +1);
The lines above save the center column of theta (which is a 2D array) into a slice of theta_middle (which is a 3D array).
Why is theta_middle a 3D array? The right hand side returns a column vector, so you should save these column vectors in a 2D array:
theta_middle(:,counter) = theta(:, N/2 +1);
The center column of array theta is also the center column of array y. Why not just save the center column of y?
Why do you save the same data the same way in 4 different arrays?
Is it guaranteed that N is even? If N is odd, you will get a non-integer array index in the lines above.
The arrays theta_middle and THETA appear to grow inside a loop. Pre-allocate them for efficiency.
The function Max_Non_linear1_IjuptilK_func() returns vector t. The t vector returned will be from the final loop pass only. This is OK, since t is the same on every pass.
The function Max_Non_linear1_IjuptilK_func() returns z, but no value is assigned to z.
  3 个评论
显示 1更早的评论
Torsten
Torsten 2022-12-2
编辑：Torsten 2022-12-2
I think you want
% theta at the middle of the layer (i.e., z=d/2)
theta_middle(:,counter) = theta(:, N/2 +1);
instead of
% theta at the middle of the layer (i.e., z=d/2)
theta_middle(counter,:,:) = theta(:, N/2 +1);
Then
theta_middle_last_time_step(counter) = theta_middle(end,counter)
University Glasgow
University Glasgow 2022-12-2
Thank you so much Torsten. This is what i want.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by