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

조회 수: 1(최근 30일)
Sebastian Sanchez 2021년 7월 27일
댓글: 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표시숨기기 이전 댓글 수: 1
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 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표시숨기기 이전 댓글 수: 1
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 2021년 7월 27일
Matlab has ode23 and ode45.
doc ode23 or ode45 for solving the differential equations.
##### 댓글 수: 0표시숨기기 이전 댓글 수: -1

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

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