Amplitude increasing in the verlet method in backward difference

Question

I am facing trouble in the increasing oscillation on a simple harmonic oscillator using a backward difference. here is my code in Scilab

function [x] = back(h, tf)
k = 2;
m = 1;
i = 2;

t(i - 1) = 0;
x(i - 1) = 10;
v(i - 1) = 0;
t(i) = t(i - 1) + h
v(i) = v(i - 1) - h * (k / m) * x(i - 1)

while t(i) < tf
    t(i + 1) = t(i) + h
    x(i + 1) = x(i - 1) - 2 * (k / m) * v(i) * h
    i = i + 1
end

plot(t, x, 'b');

endfunction

Dear Anshul, welcome to StackOverflow. it is a forum where people ask questions and others try to answer them as far as they can. It is very important that you learn how to ask good questions. You need to explain as precisely as possible that what you want to achieve, what you have tried and what your problem is. Please learn MarkDown language and code blocks using `~~~` to format your code properly. see my edits and try to follow the formating next time. — Foad S. Farimani, Mar 01 '20 at 13:06
@Anshul: your implementation of Verlet method seems wrong. Please have a look at https://en.wikipedia.org/wiki/Verlet_integration — Stéphane Mottelet, Mar 02 '20 at 14:27
@StéphaneMottelet would be nice to rectify me with the code as i am confused totally — Anshul Sharma, Mar 03 '20 at 14:49

score 2 · Accepted Answer · edited Mar 04 '20 at 16:51

2

From your code, I suppose that you are trying to implement the velocity-Verlet scheme. Here is its implementation for a simple oscillator with the differential equation:

function [x] = back(h, tf)
  k = 2;
  m = 1;

  t = 0:h:tf;
  x(1) = 10;
  v(1) = 0;

  for i=2:length(t)
    x(i) = x(i - 1) + v(i - 1) * h - k / m * x(i-1) * h^2 / 2;
    v(i) = v(i - 1) - k / m * (x(i) + x(i-1)) * h / 2;
  end

  plot(t, x, 'b');
endfunction

[x] = back(0.01, 10)

edited Mar 04 '20 at 16:51

Foad S. Farimani

12,396
15
78
193

answered Mar 04 '20 at 13:45

Stéphane Mottelet

2,853
1
9
26

1

Thank you for the edit. LaTeX support like stackexchange would be great... – Stéphane Mottelet Mar 05 '20 at 16:47
yeah sadly SO does not support LaTeX math formating. I usually use one of the LaTeX to image software like the ones [here](https://tex.stackexchange.com/a/420307/102825). [math.stackexchange](https://math.stackexchange.com/search?tab=newest&q=scilab) and [stats.stackexchange](https://stats.stackexchange.com/search?tab=newest&q=scilab) do support LaTeX and there are also som Scilab questions there which could use your kindness. Also on [Reddit](https://www.reddit.com/search/?q=scilab&sort=new) there are question. :) – Foad S. Farimani Mar 06 '20 at 13:22

score 0 · Answer 2 · answered Mar 01 '20 at 13:35

I'm not quite sure what you are trying to achieve, neither if your math is correct. But assuming that you want to solve the numerical problem of:

//coefficients of:
k = 2.;
m = 1.;

// with an initial condition of:
t(1) = 0.;
x(1) = 10.;
v(1) = 0.;

// time paramters:
N = 50;
tf = 10;
h = tf / 50.;


for ii = 2:N
    t(ii) = t(ii - 1) + h;
    x(ii) = x(ii - 1) - 2 * (k / m) * v(ii - 1) * h
    v(ii) = v(ii - 1) - h * (k / m) * x(ii - 1)
    disp(x(ii))
end


plot(t, x, 'b');

will result in:

which doesn't seem right but anyway. Please check your math again.

Amplitude increasing in the verlet method in backward difference

2 Answers2