What is wrong with my code, my equation is vdv/dy=-G*M/(y+R)^2 ?

3 次查看（过去 30 天）

显示更早的评论

Laura 2023-4-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1941144-what-is-wrong-with-my-code-my-equation-is-vdv-dy-g-m-y-r-2

编辑： Torsten 2023-4-4

cssCopy codefunction dydt = gravity(y, v, G, M, R)

% y: height above surface

% v: velocity at height y

% G: gravitational constant

% M: mass of the planet

% R: radius of the planet

dydt(1,1) = v; % v = dy/dt

dydt(2,1) = -G*M/(y+R)^2; % a = d^2y/dt^2 = dv/dt

end

scssCopy codeG = 6.67430e-11; % gravitational constant

M = 5.9722e24; % mass of Earth

R = 6.371e6; % radius of Earth

7 个评论
显示 5更早的评论隐藏 5更早的评论

Laura 2023-4-4

@Cris LaPierre Thanks, I will try now

Dyuman Joshi 2023-4-4

编辑：Dyuman Joshi 2023-4-4

@Laura I am well aware that you are trying to solve the equation by MATLAB, that is not my point.

For what time range are you trying to solve the equation? What are the initial conditions to the differential equation?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Sam Chak 2023-4-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1941144-what-is-wrong-with-my-code-my-equation-is-vdv-dy-g-m-y-r-2#answer_1208674

编辑：Sam Chak 2023-4-4

在 MATLAB Online 中打开

Hi @Laura

This is a simple example of using the solver's syntax that you can find in the ode45() documentation. If there are universal constants that they don't change, then place them inside the ode function, as suggested by @Dyuman Joshi. Otherwise, parameters that you would like to test with different value (i.e., mass), then you can pass the parameters to the ode function like @Cris LaPierre did.

What you need to do is to replace my free fall ODE with your description of Newton's law of universal gravitation.

Logically, when the object (point mass) hits the ground, then the simulation should stop. Anyhow, get familiar with the basic code. Click into the documentation for more examples.

% solving the ODE

tspan = [0 6]; % simulation time

y0 = [100 0]; % initial condition y(1) = 100 m; y(2) = 0 m/s

[t, y] = ode45(@freefall, tspan, y0);

% plotting the solution y(1)

plot(t, y(:,1), 'linewidth', 1.5),

xlabel('Time (sec)'), ylabel('Height (m)')

yline(0, '--', 'Ground level');

yregion(-80, 0)