Equation code and plots

1 次查看（过去 30 天）
Cesar Cardenas
Cesar Cardenas 2023-3-5
评论： Dyuman Joshi 2023-3-5
Hello, just wanted to know if this code would be fine for this equation? and how I could get some plots and get some error results? Any help will be greatly appreciated.thanks.
clear
clc
L = 1-0;
T = 0.06-0;
alpha = 1;
dx = 0.05;
dt = 0.001;
%Q = (alpha.*dt)/(dx.^2);
Q = 0.25;
N = L/dx +1;
M = T/dt + 1;
x = zeros(N,1);
t = zeros(M,1);
for i=1:N
x(i) = 0+(i-1).*dx;
end
for n=1:M
t(n) = 0+(n-1).*dt;
end
u = zeros(M,N);
u(:,1) = 0;
u(:,N) = 0;
u(M,2:N-1) = sin(pi.*x(2:N-1));
for n=1:M
for i=2:N-1
u(n+1,i) = Q.*u(n,i+1) + (1-2*Q).*u(n,i) + Q.*u(n,i-1)
end
end

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

回答（1 个）

Dyuman Joshi
Dyuman Joshi 2023-3-5
I did some tweaks here and there, rest of the code is good.
L = 1-0;
T = 0.06-0;
alpha = 1;
dx = 0.05;
dt = 0.001;
%Q = (alpha.*dt)/(dx.^2);
Q = 0.25;
N = L/dx + 1;
M = T/dt + 1;
%Vectorize the definition of x and t
x = (0:N-1)*dx;
t = (0:M-1)*dt;
u = zeros(M,N);
%These two lines are redundant
%u(:,1) = 0;
%u(:,N) = 0;
u(M,2:N-1) = sin(pi.*x(2:N-1));
for n=1:M
for m=2:N-1
u(n+1,m) = Q.*u(n,m+1) + (1-2*Q).*u(n,m) + Q.*u(n,m-1);
end
end
What do you want to plot? And do you have any data to find error against?
  5 个评论
显示 3更早的评论
Cesar Cardenas
Cesar Cardenas 2023-3-5
ok can you tell me what I'm missing to get the plots?
Dyuman Joshi
Dyuman Joshi 2023-3-5
Ran in:
Are use sure about this line?
u(M,2:N-1) = sin(pi.*x(2:N-1));
I think that it should be 1 instead of M, otherwise all the values in u are zero (check below)
L = 1-0;
T = 0.06-0;
alpha = 1;
dx = 0.05;
dt = 0.001;
%Q = (alpha.*dt)/(dx.^2);
Q = 0.25;
N = L/dx + 1;
M = T/dt + 1;
%Vectorize the definition of x and t
x = (0:N-1)*dx;
t = (0:M-1)*dt;
u = zeros(M,N);
u(M,2:N-1) = sin(pi.*x(2:N-1));
In the outer for loop below, if n is 1:M, then n+1 is 2:M+1, which creates an extra row in u, causing the error below in the plot. I modified it to 1:M-1.
for n=1:M-1
for m=2:N-1
u(n+1,m) = Q.*u(n,m+1) + (1-2*Q).*u(n,m) + Q.*u(n,m-1);
end
end
%Verifying that all the values in u are zero
all(u==0,'all')
ans = logical
1
figure(1)
plot (x,u(M,:)) %plot x vs u
title ('x vs u at t = 0.06')
xlabel ('x')
ylabel ('T')
figure(2)
plot(t', u(:,(N-1)/2)) %plot t vs u
%You need to take transpose of any one input
%as t is a row vector and u(:,(N-1)/2) is a column vector
title ('t vs u at x = 0.5')
xlabel ('time')
ylabel ('T')

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 2-D and 3-D Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by