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

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

[Cris LaPierre]

What error message are you seeing?

Have you copied your code correctly? These lines are not valid MATLAB systax.

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

scssCopy codeG = 6.67430e-11; % gravitational constant 
%^^^^^^^^^^^^ remove

[Laura]

We deleted the lines you said and got this error message

Error: File: untitled3.m Line: 10 Column: 1

Illegal use of reserved keyword "end".

[Dyuman Joshi]

"What is wrong with my code"

That's the thing, you have to tell us what specific problem you encountered when you ran your code.

How did you call your function? What are you attempting to solve for? v(t) or v(y)?

Note that any functions defined in a script must be at the end.

Also, define G, M and R inside the function as they are constants.

[Cris LaPierre]

I didn't say to delete them. I just think the indicated portions should not be there.

You can't have a space in a variable name.

Since it appears you are trying to solve a differential equation, have a look at this example. When you define a function, it must either be at the bottom of your script, or in its own file. See here.

So perhaps you are trying to do this? This won't do anything yet, but fixes the syntax errors I see.

G = 6.67430e-11; % gravitational constant 
M = 5.9722e24; % mass of Earth 
R = 6.371e6; % radius of Earth 

function 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  

[Laura]

@Dyuman Joshi I am trying to solve this equation using MatLab

[Laura]

@Cris LaPierre Thanks, I will try now

[Dyuman Joshi]

@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?

Answer 1

[Sam Chak]

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)

% 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  

Answer 2

[Torsten]

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])

[Dyuman Joshi]

@Torsten, the equation you wrote is incorrect.

It's

v*diff(v,y)

instead of

diff(v,y,2)

And it needs only 1 intial condition

[Torsten]

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.