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

5 visualizzazioni (ultimi 30 giorni)
Laura il 4 Apr 2023
Modificato: Torsten il 4 Apr 2023
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 CommentiMostra 6 commenti meno recentiNascondi 6 commenti meno recenti
Dyuman Joshi il 4 Apr 2023
Modificato: Dyuman Joshi il 4 Apr 2023
@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?

Accedi per commentare.

### Risposte (2)

Sam Chak il 4 Apr 2023
Modificato: Sam Chak il 4 Apr 2023
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)
% Describing the ODE
function dydt = freefall(t, y)
m = 1; % point mass
g = 9.8203; % gravity
dydt(1,1) = y(2); % dy/dt
dydt(2,1) = -g/m; % assume no air drag
end
##### 0 CommentiMostra -1 commenti meno recentiNascondi -1 commenti meno recenti

Accedi per commentare.

Torsten il 4 Apr 2023
syms G M R y v(y)
Dv = diff(v,y);
eqn = diff(v,y,2) == -G*M/(y+R)^2;
cond = [v(0)==100,Dv(0)==0];
sol = dsolve(eqn,cond)
sol =
Gnum = 6.67430e-11; % gravitational constant
Mnum = 5.9722e24; % mass of Earth
Rnum = 6.371e6; % radius of Earth
sol = vpa(subs(sol,[G M R],[Gnum,Mnum,Rnum]))
sol =
fplot(sol,[0 6])
##### 2 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Torsten il 4 Apr 2023
Modificato: Torsten il 4 Apr 2023
I also thought so first , but Laura's code (and all answers and comments thereafter) speak another language.
I think d^2v/dy^2 was misinterpreted as d(0.5*v^2)/dy.

Accedi per commentare.

### Categorie

Scopri di più su Programming in Help Center e File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by