How can I Implement Runge Kutta 4 to plot the trajectory of a given point in a vector field?
5 views (last 30 days)
Show older comments
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.
2 Comments
Accepted Answer
Chunru
on 27 Jul 2021
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
2 Comments
Chunru
on 27 Jul 2021
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.
More Answers (2)
Chunru
on 27 Jul 2021
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.
0 Comments
Chunru
on 27 Jul 2021
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
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!