Equation code and plots

15 views (last 30 days)
Cesar Cardenas
Cesar Cardenas on 5 Mar 2023
Commented: Dyuman Joshi on 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

Answers (1)

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

Sign in to comment.

Categories

Find more on Dynamic System Models in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!