RK2 Starter Scheme

조회 수: 11 (최근 30일)
Adrian Sanchez
Adrian Sanchez 2019년 12월 7일
편집: Adrian Sanchez 2019년 12월 9일
I am working on solving the following differential equation using numerical approximation using the Runge Kutta 2 Method (RK2) as a starter scheme for the leap frog method.
Consider the convection-diffusion equation in 1D:
With condtions and . Also consider the initial condition:
A. Identify if the pure convection problem gives stable solution with Forward Euler
time advancement and second-order central for spatial derivative.
B. Solve the problem for using Leapfrog for time advancement and
second-order central for spatial derivative. Use RK2 as starter scheme. Plot
the solution for t = 0,4,8. Use at least 51 points in x direction. Discuss your
solution in terms of the stability and accuracy of these schemes. Try different
appropriate values for u∆t/∆x= 0.02, 0.2, 0.8, 1.
I have the approximations and set up an anonymous function to represent the RK2 scheme and also have the leap frog method setup on paper but not yet in MATLAB.
I keep running into the error:
Index in position 1 exceeds array bounds (must not exceed51).
Error in ME_526_HW5_p1 (line 59)
k1(1,j) = dTdt(T(j-1,n), T(j,n), T(j+1,n));
I have tried stepping through the code and I'm not sure where I am going wrong.
I have also tried the initial conditions as follows:
% Initialization of T, k1, and k2
T = zeros(length(x),length(t));
k1 = zeros(length(x),length(t));
k2 = zeros(length(x),length(t));
And recieve a similar error.
I have also tried the for loop as follows:
for n = 1:length(t)
if n == 1 % Calculate intial condition
j = 1;
while x(j) <= 0.2
T(j,n) = 1-(10*x(j)-1)^2;
j = j+1;
end
elseif n == 2 % RK2 Starter Scheme
for j = 1:length(x)
k1(1,j) = dTdt(T(j-1,n), T(j,n), T(j+1,n));
end
% code that will be polulated by k2 and calculate the time step we are interseted in
else % Will have the bulk of the problem (leap frog method)
end
end
But then an index of zero occurs and MATLAB can't handle that.
The relavent code is provided below.
u = 0.08
alpha = 0
dx = .02;
dt = .005;
x = 0:dx:1;
t = 0:dt:8;
% Anonymous Function
dTdt = @(T0, T1, T2) -u/(2*dx)*(T2 - T0) + alpha/(2*dx^2)*(T2 - 2*T1 + T0);
% Initialization of T, k1, and k2
T = zeros(length(x),length(t));
k1 = zeros(1,length(x));
k2 = zeros(1,length(x));
% Intital Condition
% T(x,0) = 1-(10*x-1)^2 for 0 <= x <= 0.2
% T(1,t) = 0
% T(x,0) = 0
for n = 1:length(t)
if n == 1 % Calculate intial condition
j = 1;
while x(j) <= 0.2
T(j,n) = 1-(10*x(j)-1)^2;
j = j+1;
end
elseif n == 2
for j = 2:length(x)
k1(1,j) = dTdt(T(j-1,n), T(j,n), T(j+1,n));
end
% code that will be polulated by k2 and calculate the time step we are interseted in
else % Will have the bulk of the problem (leap frog method)
end
end
  댓글 수: 2
J. Alex Lee
J. Alex Lee 2019년 12월 8일
There's not really a question in there, but to address the seemingly most salient aspect of your issue, the zero-index, you seem to be attempting to apply a central difference to a boundary node (in particular the left boundary) so that you are trying to access the "first minus one" node, which doesn't exist (1-1=0).
It looks to me that your first task should be to undertand how to deal with boundary conditions in a finite difference scheme - do you have a handle on how to solve ODEs (diffeq's in 1 variable, such as steady state version of your problem) as BVPs (high order ODEs with boundary conditions on both sides), nevermind PDEs (diffeq's in multiple variables, as you have here)?
Adrian Sanchez
Adrian Sanchez 2019년 12월 9일
편집: Adrian Sanchez 2019년 12월 9일
Thank you for your response, though I did end up figuring out what the issue was. It was that I was ending my loop for k1 at j = 51 which could not be done because I was also referencing j = 52 in my loop which exceeded the limits of the matrix T. The loops to solve the RK2 starter scheme should have taken the following form:
for j = 2:length(x)-1
k1(j-1) = dTdt(T(j-1,n-1), T(j,n-1), T(j+1,n-1));
end
for j = 2:length(x)-1
k2(j)= dTdt(dt(i)*(T(j-1,n-1) + k1(j-1)/2), dt(i)*(T(j,n-1) + ...
k1(j)/2), dt(i)*(T(j+1,n-1)+ k1(j+1)/2));
T(j,n) = T(j,n-1) + k2(j);
end
I do understand how to solve simple ODE's, my issue was a MATLAB syntax that I was failing to recognize.

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

답변 (0개)

카테고리

Help CenterFile Exchange에서 Particle & Nuclear Physics에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by