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일

0 개 추천

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일
Thank you, I would like to plot u, but not sure how to do ti? no, I do not have any data.
This is what I'm trying to plot, but I'm getting errors.
Error using plot
Vectors must be the same length.
Error in HW1 (line 31)
plot(t, u(:,(N-1)/2)) %plot t vs u
As always any help wil be appreciated.
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
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
title ('t vs u at x = 0.5')
xlabel ('time')
ylabel ('T')
Dyuman Joshi
Dyuman Joshi 2023년 3월 5일
plot U with respect to what? Generally, you need x and y values for a plot.
If you don't have any data, how do you want to calculate the error?
Cesar Cardenas
Cesar Cardenas 2023년 3월 5일
ok can you tell me what I'm missing to get the plots?
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')

댓글을 달려면 로그인하십시오.

카테고리

도움말 센터File Exchange에서 Dynamic System Models에 대해 자세히 알아보기

태그

질문:

2023년 3월 5일

댓글:

2023년 3월 5일

Community Treasure Hunt

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

Start Hunting!

Translated by