Runge-Kutta integration of the Taylor-Maccoll eq

Hi, I'm trying to integrate the Taylor-Maccoll equation in order to determine the velocity flowfield in a conical shockwave, moving at hypersonic speed. As explained in the photo ive attached, the equation can be represented as a standard ODE with solution vector. I have initial conditions for both vr and vr', which I am then supposed to place into the second vector, y'. I then integrate y' using the Runge-Kutta method and at some point in the integration vr' should equal (roughly) 0. I have a step size of h, which defines the increment in theta (polar coordinates) I take from the shockwave, towards the cone angle.

I am aware of two codes which have been uploaded to the MATLAB servers which I can download, however these codes do not include the discrete values of velocity at every angle of theta, and so I must code these a different way.

In my code I have manually typed up the Runge-Kutta method but I think there is some problem in my code as my vrdash values are not equating to zero at the correct theta value (noweher near in fact!).

% initial values for velocity just after the shock
vr = 0.9078; 
vrdash = -0.1262;
theta_s = 14.35;

%initial y matrix
y1 = zeros(1, length(x));
y2 = zeros(1, length(x));

%initial conditions for velocity
y1(1) = vr;
y2(1) = vrdash;

%defining the interval and step size (thetas_s = half angle of shock wave)
h = -0.01;           %negative step size as moving towards zero in theta direction
x = theta_s:h:0;     %integrate between shock angle and zero

g = 1.4;             %ratio of specific heats
       
t(1) = 0;

%defining function from y' equation in photo
A = (g-1)/2;  
    F1 = @(t, y1, y2) y2;
    F2 = @(t, y1, y2) (y1*y2^2-A*(1-y1^2-y2^2)*(2*y1+y2*cotd(theta_s)))...
                      /(A*(1-y1^2-y2^2) - y2^2);

%initial loop for Runge-Kutta method
 for i=1:length(x)  
   
    k1 = h*F1(t, y1(i), y2(i));
    m1 = h*F2(t, y1(i), y2(i));
    
    k2 = h*F1(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
    m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*m1);
    
    k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
    m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*m2));
    
    k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
    m4 = h*F2(t+h, (y1(i)+h),(y2(i)+m3*h));
 
    % main equation
    y1(i+1) = y1(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;
    y2(i+1) = y2(i) + (1/6)*(m1+2*m2+2*m3+m4)*h;             
 end  
 end

I'm sorry if I havent explained this properly, if there's anything i'm missing I'd be happy to explain. Any help would be greatly appreciated as I've been stuck on this code for nearly a month now. Thanks!

Answer 1

[Alan Stevens]

Shouldn't

   
    k1 = h*F1(t, y1(i), y2(i));
    m1 = h*F1(t, y1(i), y2(i));
    
    k2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
    m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
    
    k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
    m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
    
    k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
    m4 = h*F2(t+h, (y1(i)+h),(y2(i)+k3*h));
 

be

    k1 = h*F1(t, y1(i), y2(i));
    m1 = h*F2(t, y1(i), y2(i));
    
    k2 = h*F1(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
    m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*m1);
    
    k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
    m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*m2));
    
    k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
    m4 = h*F2(t+h, (y1(i)+h),(y2(i)+m3*h));

[Robert Wake]

Hi Alan, yes i think you're right. I've copied this in, however the code is still not working in the desired way

[Alan Stevens]

LIke so (you had too many h's!).

% initial values for velocity just after the shock
vr = 0.9078; 
vrdash = -0.1262;
theta_s = 14.35;

%defining the interval and step size (thetas_s = half angle of shock wave)
h = -0.01;           %negative step size as moving towards zero in theta direction
x = theta_s:h:0;     %integrate between shock angle and zero
                %%%%%%  DEFINE x BEFORE USING IT TO CONSTRUCT y

%initial y matrix  
y1 = zeros(1, length(x));
y2 = zeros(1, length(x));

%initial conditions for velocity
y1(1) = vr;
y2(1) = vrdash;

g = 1.4;             %ratio of specific heats
       
t(1) = 0;         %%%%% t IS NEVER USED

%defining function from y' equation in photo
A = (g-1)/2;  
    F1 = @(t, y1, y2) y2;
    F2 = @(t, y1, y2) (y1*y2^2-A*(1-y1^2-y2^2)*(2*y1+y2*cotd(theta_s)))...
                      /(A*(1-y1^2-y2^2) - y2^2);

%initial loop for Runge-Kutta method
 for i=1:length(x)-1  
   
    k1 = h*F1(t, y1(i), y2(i));
    m1 = h*F2(t, y1(i), y2(i));
    
    k2 = h*F1(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
    m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*m1);
    
    k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
    m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*m2));
    
    k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
    m4 = h*F2(t+h, (y1(i)+h),(y2(i)+m3*h));
 
    % main equation   NO NEED FOR ANY MORE h's
    y1(i+1) = y1(i) + (1/6)*(k1+2*k2+2*k3+k4);
    y2(i+1) = y2(i) + (1/6)*(m1+2*m2+2*m3+m4);             
 end  
 
 plot(x,y1),grid
 xlabel('x'),ylabel('y1')

[Robert Wake]

Alan thank you so much! really appreciate it!