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

Question

Sebastian Sanchez 2021-7-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/886549-how-can-i-implement-runge-kutta-4-to-plot-the-trajectory-of-a-given-point-in-a-vector-field

评论： Chunru 2021-7-27

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 个评论
显示无隐藏无

James Tursa 2021-7-27

What have you tried so far?

Sebastian Sanchez 2021-7-27

在 MATLAB Online 中打开

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 2021-7-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/886549-how-can-i-implement-runge-kutta-4-to-plot-the-trajectory-of-a-given-point-in-a-vector-field#answer_754419

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Sebastian Sanchez 2021-7-27

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 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.

请先登录，再进行评论。