How can I Implement Runge Kutta 4 to plot the trajectory of a given point in a vector field?

Hello there,

I'm trying to approximate the trajectory of a particle with initial coordinates a,b and an initial velocity vector <P,Q> on a vector field of the form f(x,y) = <u(x,y), v(x,y)>.

Is there any form to calculate this using RK4? The output should be a plot of the vector field and the trajectory of the particle.

I have tried everything but I still got no solution to this problem, maybe you guys could tell me a different approach to solve this using Matlab.

Thanks for reading at this point and attending my request, have a nice day.

[James Tursa]

What have you tried so far?

[Sebastian Sanchez]

I have a vector field in this form

clc;
clear;
close all;
[a,b] = meshgrid(0:0.2:2,0:0.2:2);

u = cos(a).*b;
v = sin(a).*b;


izq=0;
der=2;
h=0.2;
t=(izq:h:der);
x=zeros(length(t),1);
y=zeros(length(t),1);
syms l m;
f=@(t,x);
x(1)=0;
for i=1:length(t)-1
    for j=1:length (t)-1
        k1=f(t(i),x(i));
        k2=f(t(i)+h/2,x(i)+h/2*k1);
        k3=f(t(i)+h/2,x(i)+h/2*k2);
        k4=f(t(i)+h,x(i)+h*k3);
        y(i+1)=y(i)+(1/6)*h*(k1+2*k2+2*k3+k4);
    end
end

figure
plot(t,y);
hold on
quiver(a,b,u,v);
hold off

Im looking for a differential equation (f in the code)that goes along the direction of the vectors

Answer 1

[Chunru]

As you have the intial point (0,0) and the intial speed (0,0), so the solution never progress to a different point. If you change the speed formla or initial point, then you can see something different.

[a,b] = meshgrid(0:0.2:2,0:0.2:2);
u = cos(a).*b;
v = sin(a).*b;

izq=0;
der=2;
h=0.2;
t=(izq:h:der);
%x0 = [0; 0];
x0 = [0.2; 0.2];

[t, y] = ode23(@odefun,t,x0);
figure
%plot(t,y);
plot(y(:,1), y(:,2))
hold on
quiver(a,b,u,v);
hold off

function y = odefun(t, x)
%u = cos(a).*b;
%v = sin(a).*b;
y(1,1) = cos(x(1)) .* x(2);       % u(x, y)
y(2,1) = sin(x(1)) .* x(2);       % v(x, y)
end

[Sebastian Sanchez]

This solution is the one i need. Although, im unsure of how the part [t, y] = ode23(@odefun,t,x0) works, how ode23 gets the parameters that should be passed in odefun (t and x)?

[Chunru]

Read https://blogs.mathworks.com/cleve/2014/05/26/ordinary-differential-equation-solvers-ode23-and-ode45/

For an ODE , you need to specify the function f(t, x). For 2D ODE, you specify f(t, x, y). In [t, y] = ode23(@odefun,t,x0) , @(odefun) define f(t, x, y) as a function handle; t is the time points and x0 is the initial value. Then ode23 is doing something similar to Runge Kutta you have implemented.

Answer 2

[Chunru]

Matlab has ode23 and ode45.

https://blogs.mathworks.com/cleve/2014/05/26/ordinary-differential-equation-solvers-ode23-and-ode45/

doc ode23 or ode45 for solving the differential equations.

Answer 3

[Chunru]

tspan = 0:.1:10;
x0 = [0; 0];
[t,x] = ode23(@odefun,tspan,x0);
subplot(211); plot(t, x); xlabel('t'); legend('x', 'y')
subplot(212); plot(x(:,1), x(:,2)); xlabel('x'); ylabel('y')

function y = odefun(t, x)
y(1,1) = (x(1).^2 + x(2).^2) * cos(x(1)) + .1;       % u(x, y)
y(2,1) = (x(1).^2 + x(2).^2) * sin(x(2)) + .2;       % v(x, y)
end