Equation code and plots

1 vue (au cours des 30 derniers jours)
Cesar Cardenas
Cesar Cardenas le 5 Mar 2023
Commenté : Dyuman Joshi le 5 Mar 2023
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

Réponses (1)

Dyuman Joshi
Dyuman Joshi le 5 Mar 2023
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 commentaires
Cesar Cardenas
Cesar Cardenas le 5 Mar 2023
ok can you tell me what I'm missing to get the plots?
Dyuman Joshi
Dyuman Joshi le 5 Mar 2023
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')

Connectez-vous pour commenter.

Catégories

En savoir plus sur Dynamic System Models 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!

Translated by