필터 지우기
필터 지우기

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
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일
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 CenterFile 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