I am trying to use a Runga Kutta 4th order method to solve for the Taylor Maccoll equation.

I have been given instructions to make this code and I am stuck on the following,

" For this section we are interested in when the values change from negative to positive so we will be using a “while” loop. This will keep performing the loop continuously until the defined condition is met. Because we are interested in when v_theta becomes positive we will have our while loop continue as long as v_theta is <0. We can see on the bottom of the given page on Runge-Kutta that w2 = u2 so the start of our loop will be “while w2<0”. Once we begin our loop we are essentially copy pasting the exact equations from the page on Runge-Kutta directly into our code. And then before we end our loop we have to update our old variables to make sure we keep iterating . So theta has to be increased by h (delta theta) and i has to be increased by 1."

These are the equations that were given:

And this is the code I have so far:

clc; clear all;
%Ts-shock angle, Tc-Cone angle assumption, Td- delta,
Ts=17;
Tc=0;
M1=5;
k=1.4;
Td=atand(2*cotd(Ts)*(((M1^2)*((sind(Ts))^2)-1)/((M1^2)*(k+cosd(2*Ts))+2)));
Mn1=M1*sind(Ts);
N=Mn1^2+(2/(k-1));
D=((2*k)/(k-1))*(Mn1^2)-1;
Mn2=sqrt(N/D);
M2=Mn2/(sind(Ts-Td));
%V-total vel, Vr- radial vel, VT-theta vel
V=sqrt(((k-1)/2)*(M2^2)/(1+((k-1)/2)*(M2^2)));
VT=-V*sind(Ts-Td);
Vr=V*cosd(Ts-Td);

h=-0.1;                                                                   
T=Ts:h:Tc;                                                          
A=(k-1)/2;
u1=zeros(1,length(T));
u2=zeros(1,length(T));
u1(1)=Vr;
u2(1)=VT;
T(1)=Ts; 
F1=@(T, u1, u2) u2;
F2=@(T, u1, u2) (u1*(u2^2)-A*(1-(u1^2)-(u2^2))*(2*u1+u2*cotd(Ts)))/(A*(1-(u1^2)-(u2^2))-(u2^2));
for i=1:(length(T)-1)
    while u2<0
    k11=h*F1(T, u1(i), u2(i));
    k12=h*F2(T, u1(i), u2(i));    
    k21=h*F1(T+0.5*h, u1(i)+0.5*k11, u2(i)+0.5*k12);
    k22=h*F2(T+0.5*h, u1(i)+0.5*k11, u2(i)+0.5*k12);    
    k31=h*F1(T+0.5*h, u1(i)+0.5*k21, u2(i)+0.5*k22);
    k32=h*F2(T+0.5*h, u1(i)+0.5*k21, u2(i)+0.5*k22);    
    k41=h*F1(T+h, u1(i)+k31,u2(i)+k32);
    k42=h*F2(T+h, u1(i)+k31,u2(i)+k32);
    u1(i+1)=u1(i)+(1/6)*(k11+2*k21+2*k31+k41);
    u2(i+1)=u2(i)+(1/6)*(k12+2*k22+2*k32+k42);
    end
end

When I run the code the first value of u2 is negative but immediately after it becomes positive. I am not sure if I have written the runga kutta method wrong or if the initial conditions are wrong. My other suspicion is that on the line where I define F2, Ts should be T, but I get dimension errors and I don't know if thats the right way to go and how to fix it.

[Dyuman Joshi]

"When I run the code the first value of u2 is negative"

Because you have defined it as negative.

VT=-V*sind(Ts-Td);
%  ^ minus sign right here
u1(1)=Vr;

I don't understand why you have used a while loop?

Yes, in the definition of F2, it should be T. And you need to use the ith element of T in each corresponding step.

clc; clear all;
%Ts-shock angle, Tc-Cone angle assumption, Td- delta,
Ts=17;
Tc=0;
M1=5;
k=1.4;
Td=atand(2*cotd(Ts)*(((M1^2)*((sind(Ts))^2)-1)/((M1^2)*(k+cosd(2*Ts))+2)));
Mn1=M1*sind(Ts);
N=Mn1^2+(2/(k-1));
D=((2*k)/(k-1))*(Mn1^2)-1;
Mn2=sqrt(N/D);
M2=Mn2/(sind(Ts-Td));
%V-total vel, Vr- radial vel, VT-theta vel
V=sqrt(((k-1)/2)*(M2^2)/(1+((k-1)/2)*(M2^2)));
VT=-V*sind(Ts-Td);
Vr=V*cosd(Ts-Td);d

h=-0.1;                                                                   
T=Ts:h:Tc;                                                          
A=(k-1)/2;
u1=zeros(1,length(T));
u2=zeros(1,length(T));
u1(1)=Vr;
u2(1)=VT;
T(1)=Ts; 
F1=@(T, u1, u2) u2;
F2=@(T, u1, u2) (u1*(u2^2)-A*(1-(u1^2)-(u2^2))*(2*u1+u2*cotd(T)))/(A*(1-(u1^2)-(u2^2))-(u2^2));
for i=1:(length(T)-1)
    %while u2(i)<0
    k11=h*F1(T(i), u1(i), u2(i));
    k12=h*F2(T(i), u1(i), u2(i));    
    k21=h*F1(T(i)+0.5*h, u1(i)+0.5*k11, u2(i)+0.5*k12);
    k22=h*F2(T(i)+0.5*h, u1(i)+0.5*k11, u2(i)+0.5*k12);    
    k31=h*F1(T(i)+0.5*h, u1(i)+0.5*k21, u2(i)+0.5*k22);
    k32=h*F2(T(i)+0.5*h, u1(i)+0.5*k21, u2(i)+0.5*k22);    
    k41=h*F1(T(i)+h, u1(i)+k31,u2(i)+k32);
    k42=h*F2(T(i)+h, u1(i)+k31,u2(i)+k32);
    u1(i+1)=u1(i)+(1/6)*(k11+2*k21+2*k31+k41);
    u2(i+1)=u2(i)+(1/6)*(k12+2*k22+2*k32+k42);
    %end
end
u1

u1 = 1×171

0.8730 0.8797 0.8649 0.8283 0.7757 0.7170 0.6527 0.5838 0.5109 0.4352 0.3576 0.2794 0.2017 0.1263 0.0550 -0.0093 -0.0583 -0.0815 -0.0963 -0.1165 -0.1501 -0.2065 -0.2596 -0.3033 -0.3167 -0.3196 -0.3160 -0.3025 -0.2513 -0.1712

u2

u2 = 1×171

-0.1486 0.0166 0.2570 0.4928 0.5575 0.6160 0.6673 0.7105 0.7445 0.7684 0.7813 0.7816 0.7679 0.7372 0.6842 0.5929 0.4422 0.1280 0.1703 0.2429 0.3761 0.6082 0.3552 0.4034 0.0533 0.0005 -0.0778 -0.2045 -0.8070 -0.7913

[Torsten]

Your definition of F2 must read

F2=@(T, u1, u2) (u1*(u2^2)-A*(1-(u1^2)-(u2^2))*(2*u1+u2*cotd(T)))/(A*(1-(u1^2)-(u2^2))-(u2^2));

instead of

F2=@(T, u1, u2) (u1*(u2^2)-A*(1-(u1^2)-(u2^2))*(2*u1+u2*cotd(Ts)))/(A*(1-(u1^2)-(u2^2))-(u2^2));

And your independent variable T must be set up in radians, not in degrees. Thus cotd must be changed to cot.

And be careful: You integrate backwards. For this, you have to make a change of variables

y(t) = ytilde(ttilde), ttilde = -t

and your equation

dy/dt = f(t,y) , y(t0) = y0

changes to

dytilde/dttilde = -f(-ttilde,ytilde(ttilde)), ytilde(-t0) = y0

[Jiho Han]

u2 is negative until it reaches the the cone surface which is why it should be negative until T(theta) is around 11.5 degrees. Which is why I am confused to why it immediately becomes positive.

[Torsten]

So you already corrected all the errors I listed ?

Answer 1

[Nipun]

Hi Jiho,

I understand that you are trying to solve the system of equations using fourth order Range Kutta method in MATLAB.

Based on the provided information and the attached code, I suggest the following amendments to your code:

You've correctly identified that Ts should be T(i) in the definition of F2, as you're dealing with a changing angle over the iterations. Moreover, the while loop should be inside the for loop, and you should be checking the condition for u2(i) instead of u2, which would otherwise reference the entire array.

Additionally, the while loop should not be running for the entire length of T. Instead, it should only run until the condition u2(i) < 0 is met. Once u2(i) becomes positive, the while loop should stop, and the for loop should continue to the next iteration.

Here's a revised version of the loop structure with the necessary changes:

for i = 1:(length(T)-1)
    % Runge-Kutta calculations for the i-th step
    % The while loop is used to continue calculations until u2(i) becomes positive
    while u2(i) < 0
        k11 = h*F1(T(i), u1(i), u2(i));
        k12 = h*F2(T(i), u1(i), u2(i));    
        k21 = h*F1(T(i)+0.5*h, u1(i)+0.5*k11, u2(i)+0.5*k12);
        k22 = h*F2(T(i)+0.5*h, u1(i)+0.5*k11, u2(i)+0.5*k12);    
        k31 = h*F1(T(i)+0.5*h, u1(i)+0.5*k21, u2(i)+0.5*k22);
        k32 = h*F2(T(i)+0.5*h, u1(i)+0.5*k21, u2(i)+0.5*k22);    
        k41 = h*F1(T(i)+h, u1(i)+k31, u2(i)+k32);
        k42 = h*F2(T(i)+h, u1(i)+k31, u2(i)+k32);
        
        u1(i+1) = u1(i) + (1/6)*(k11 + 2*k21 + 2*k31 + k41);
        u2(i+1) = u2(i) + (1/6)*(k12 + 2*k22 + 2*k32 + k42);
        
        % Update the angle and the index for the next iteration
        T(i+1) = T(i) + h;
        
        % Break the while loop if u2(i+1) is non-negative
        if u2(i+1) >= 0
            break;
        end
        
        i = i + 1; % Increment the index
        if i >= length(T)
            break; % Avoid exceeding the bounds of the arrays
        end
    end
    
    % Check if the while loop was broken due to u2(i+1) being non-negative
    if u2(i+1) >= 0
        break; % Exit the for loop early since the condition is met
    end
end

Make sure to also correct the definition of F2 to use T(i) instead of Ts. Here's the corrected line:

F2 = @(T, u1, u2) (u1*(u2^2)-A*(1-(u1^2)-(u2^2))*(2*u1+u2*cotd(T)))/(A*(1-(u1^2)-(u2^2))-(u2^2));

Remember to update the function handle F2 wherever it is used in the Runge-Kutta calculations to pass T(i) instead of T.

Here is a similar query on implementing Range Kutta method in MATLAB: Runge-Kutta 4th order method - MATLAB Answers - MATLAB Central (mathworks.com)

Hope this helps.

Regards,

Nipun