I am trying to use a Runga Kutta 4th order method to solve for the Taylor Maccoll equation.
10 views (last 30 days)
Show older comments
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1359348/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1359353/image.png)
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.
4 Comments
Answers (1)
Nipun
on 29 Dec 2023
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
0 Comments
See Also
Categories
Find more on Calculus 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!