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

조회 수: 28 (최근 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
Sebastian Sanchez
Sebastian Sanchez 2021년 7월 27일
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

댓글을 달려면 로그인하십시오.

채택된 답변

Chunru
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
Sebastian Sanchez
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
Chunru 2021년 7월 27일
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
Chunru 2021년 7월 27일
Matlab has ode23 and ode45.
doc ode23 or ode45 for solving the differential equations.

Chunru
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

제품


릴리스

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by