How can I Implement Runge Kutta 4 to plot the trajectory of a given point in a vector field?
5 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Chunru
2021-7-27
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 个评论
Chunru
2021-7-27
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.
更多回答(2 个)
Chunru
2021-7-27
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 个评论
Chunru
2021-7-27
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 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!