I am trying to create a Dufort Frankel Finite Difference Scheme. I keep getting an error that says "Index in position 2 is invalid. Array indices must be positive integers or logical values." I think I understand that Matlab must start at 1 and that my tdx-1 messes this up. I have tried several things to change this but I am having trouble getting it to work. Could someone inform me of how to correct this issue?
%%%%% dT/dt = d^2T/dx^2 - partial derivatives
%%%% T = temperature = f(x,t)
%%%%% time derivative = second spatial derivative
%%%%(T(x,t+dt) - T(x,t))/dt = (T(x+dx,t) - 2*T(x,t) + T(x-1,t))/ dx^2
%%% Propagate this system in time
%%%% Solve for T(x,t+dt) =T(x,t) + k*dt/dx^2 * (T(x+dx,t) - 2*T(x,t) +
%%%% T(x-1,t))
%%%% length of the pipe
k = 2;
L = 10;
N = 10;
x_vec = linspace(0,L,N);
dx = x_vec(2)-x_vec(1);
dt = 0.5*(dx^2)/(2*k);
t_vec = 0:dt:10;
T_mat = zeros(length(x_vec),length(t_vec));
T_mat(1,:) = 200;
T_mat(end,:) = 150;
for tdx = 1;length(t_vec)-1;
for idx = 2:length(x_vec)-1
T_mat(idx,tdx+1) = T_mat(idx,tdx-1) + 2*k*dt/dx^2*( T_mat(idx+1,tdx) - (T_mat(idx,tdx+1)+T_mat(idx,tdx-1)) + (T_mat(idx-1,tdx)));
end
end
[tt,xx] = meshgrid(t_vec,x_vec);
mesh(xx,tt,T_mat)
xlabel('X coordinate (m)')
ylabel('Time (sec)')
zlabel('Temperature (F)')

 採用された回答

Cris LaPierre
Cris LaPierre 2019 年 3 月 11 日
編集済み: Cris LaPierre 2019 年 3 月 12 日

0 投票

The ";" instead of ":" is one issue as madhan ravi pointed out. However, another issue is with your indexing of tdx in your for loop.
for tdx = 1:length(t_vec)-1
for idx = 2:length(x_vec)-1
T_mat(idx,tdx+1) = T_mat(idx,tdx-1) + 2*k*dt/dx^2*( T_mat(idx+1,tdx) - (T_mat(idx,tdx+1)+T_mat(idx,tdx-1)) + (T_mat(idx-1,tdx)));
end
end
Note that your 2nd index for T_mat is sometimes tdx-1. But you start with tdx=1. As you point out, MATLAB indexes starting at 1, so this is throwing an error.
Did you perhaps get your idx and tdx backwards? You start with idx=2.

13 件のコメント

Hunter Sylvester
Hunter Sylvester 2019 年 3 月 12 日
I see what you are saying. I fixed the colon and semicolon mistake then I changed tdx = 2 and idx = 2. This has actually produced something, thank you very much. The graph produced is not very good how can I correct this or make it look more denounced to the viewer.
This is what you meant to do right?
for tdx = 2:length(t_vec)-1
for idx = 2:length(x_vec)-1
T_mat(idx,tdx+1) = T_mat(idx,tdx-1) + 2*k*dt/dx^2*( T_mat(idx+1,tdx) - (T_mat(idx,tdx+1)+T_mat(idx,tdx-1)) + (T_mat(idx-1,tdx)));
end
end
Hunter Sylvester
Hunter Sylvester 2019 年 3 月 12 日
The Dufort Frankel scheme is the image attached below. What should I change the k and dt to equal so that this works properly.
madhan ravi
madhan ravi 2019 年 3 月 12 日
編集済み: madhan ravi 2019 年 3 月 12 日
I have no idea how this was different from the answer below.
Cris LaPierre
Cris LaPierre 2019 年 3 月 12 日
I'm not sure what this is supposed to look like. You know more about it than I do.
To explore a parameter space, I'd suggest putting your code into a live script and adding slider controls for k and dt. That way you can adjust the values until you get what you're expecting, and can then see what values work.
finiteDiffControls.png
Hunter Sylvester
Hunter Sylvester 2019 年 3 月 12 日
Wow, did not realize you could do that. Thank you so much.
Cris LaPierre
Cris LaPierre 2019 年 3 月 12 日
For madhan ravi - As you pointed out, the semicolon was wrong, but not the source of the error message. It was necessary to fix the indexing of T_mat for when tdx=1:
T_mat(idx,tdx-1)
Hunter Sylvester
Hunter Sylvester 2019 年 3 月 12 日
I am looking for the insert option along with the slider option in the insert tab. Would I have this option if I have the student version of Matlab? I only bought the $99 version.
Cris LaPierre
Cris LaPierre 2019 年 3 月 12 日
編集済み: Cris LaPierre 2019 年 3 月 12 日
There is no difference is capabilities in Student version. It would just depend on what release you purchased. Interactive controls were introduced in r2018a I believe.
madhan ravi
madhan ravi 2019 年 3 月 12 日
編集済み: madhan ravi 2019 年 3 月 12 日
@Cris: So that’s what I mentioned that tdx iterates from 2?
Cris LaPierre
Cris LaPierre 2019 年 3 月 12 日
Ah, got it. FWIW, we were writing our responses at the same time. I guess the difference here then was that this answer explained the error, a helpful addition for new users who don't understand why they are getting an error.
madhan ravi
madhan ravi 2019 年 3 月 12 日
編集済み: madhan ravi 2019 年 3 月 12 日
Hmm yes that‘a true , by the way the slider control is amazing, didn’t know this feature.Any links to the documentation would be great.
Cris LaPierre
Cris LaPierre 2019 年 3 月 12 日
I put this link above. Just sliders and drop downs in 18a and 18b, but more coming soon!
madhan ravi
madhan ravi 2019 年 3 月 12 日
Cool!

サインインしてコメントする。

その他の回答 (1 件)

madhan ravi
madhan ravi 2019 年 3 月 11 日

0 投票

for tdx = 2:length(t_vec)-1

カテゴリ

ヘルプ センター および File ExchangeProgramming についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by