Thread: Fourth Order Runge-Kutta for two-simultaneous 1st ODE

I am having real trouble in creating a fourth order Runge-Kutta method for a system of two simultaneous first order ordinary differential equations.

It is supposed to model the wake off an aircraft wing, in terms of it's x and y coordinates.


The wing has been split into 10 parts, i.e. 5 on each side of the wing. And the program should give 10 values of x and y for each time step (i.e. 10 lots of 500 time steps in this case)
Any help with this problem would be gratefully appreciated..My attempted matlab program is below, but i fear ive got lost....

clc
m=10;
m2=m/2;
dx=1/m;
NMAX=500;
nfreq=1;
t=0;
x(1)=-1+dx;
ndt=0;
h=1;
dt=1/25/m;

for J=2:m;
x(J)=x(J-1)+2*dx;


end
for J=1:m2;
y(J)=0;
g(J)=((1-((x(J)+dx)^2))^0.5)-((1-((x(J)-dx)^2))^0.5);

end

m2p1=m2+1;

for J=m2p1:m;
y(J)=0;
g(J)=-g(m+1-J);

end


U=0;
V=0;
XI=-0.9;
YI=0;

NDT=0;
T=0;

while NMAX>NDT
NDT=NDT+1;
T=T+dt;

for J=1:10;
if J~=I;

U1=U+g(J)*(YI-Y(J))/((XI-X(J))^2 + (YI-Y(J))^2);
V1=V-g(J)*(XI-X(J))/((XI-X(J))^2 + (YI-Y(J))^2);
D1X= dt*U1;
D1Y= dt*V1;
X2=XI+D1X/2;
Y2=YI+D1Y/2;
U2=U+g(J)*(Y2-Y(J))/((X2-X(J))^2 + (Y2-Y(J))^2);
V2=V-g(J)*(X2-X(J))/((X2-X(J))^2 + (Y2-Y(J))^2);
D2X= dt*U2;
D2Y= dt*V2;
X3=XI+D2X/2;
Y3=YI+D2Y/2;
U3=U+g(J)*(Y3-Y2)/((X3-X2)^2 + (Y3-Y2)^2);
V3=V-g(J)*(X3-X2)/((X3-X2)^2 + (Y3-Y2)^2);
D3X= dt*U3;
D3Y= dt*V3;
X4=XI+D3X;
Y4=YI+D3Y;
U4=U+g(J)*(Y4-Y3)/((X4-X3)^2 + (Y4-Y3)^2);
V4=V-g(J)*(X4-X3)/((X4-X3)^2 + (Y4-Y3)^2);
D4X= dt*U4;
D4Y= dt*V4;

X(J+1) = X(J)+(D1X+2*D2X+2*D3X+D4X)/6.0;
Y(J+1) = Y(J)+(D1Y+2*D2Y+2*D3Y+D4Y)/6.0;

X(I)=X(J+1);
Y(I)=Y(J+1);
end
end
for I=m2p1:m;
X(I)=-X(m+1-I);
Y(I)=-Y(m+1-I);

end

X(1:10)
Y(1:10)

end

g(1:10)

Originally Posted by p123nky

I am having real trouble in creating a fourth order Runge-Kutta method for a system of two simultaneous first order ordinary differential equations.

It is supposed to model the wake off an aircraft wing, in terms of it's x and y coordinates.


The wing has been split into 10 parts, i.e. 5 on each side of the wing. And the program should give 10 values of x and y for each time step (i.e. 10 lots of 500 time steps in this case)
Any help with this problem would be gratefully appreciated..My attempted matlab program is below, but i fear ive got lost....

Introduce the state vector .

Now write a function that returns .

You should now be able to use the update equations for the scalar RK4 on these vectors without modification.

CB