Simple harmonic motion, RK4 method, matlab plot












0












$begingroup$


So I am having a issue plotting a simply harmonic motion of the form



$$frac{d^2y}{dx^2}+frac{k}{m}y=0$$



Using the RK4 method in matlab. The issue is that my code is not producing the expected plotted and I am not entirely sure, if it my RK4 that is wrong or my actual code that is wrong.



So the solution to the above equation is $y=3cos(frac{k}{m}x)$, if i set k=1 and m=1, I should produce the following graph.



enter image description here



Using Demos



The graph that is being produced by my code is



enter image description here



here is the code which produces the given graph



clc;
clear all;
%Defing intial paramters
%dt=0.001;

h=0.001;
k=1;
m=1;
t=0:h:100;

x=zeros(length(t),1);
v=zeros(length(t),1);
x(1)=3;
v(1)=0;

%defing functions

f=@(t,x,v) v;
g=@(t,x,v) -k/m*x;





for n=1:(length(t)-1);
k1x=f(t(n),x(n),v(n));
k1v=g(t(n),x(n),v(n));

k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k3x=f(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);
k3v=g(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);

x(n+1)=x(n)+h/6*(k1x+2*k2x+2*k3x+k4x);
v(n+1)=v(n)+h/6*(k1v+2*k2v+2*k3v+k4v);


end

plot(t,x);


`
I cannot understand to why my graph is being produced this way, I have gone through the code a couple of time, but cant see to find what is wrong. My only conclusion so far is that it either a variable I have missed or place wrong or my for stamen t is incorrect but cant see that being the case due to, using the equation for RK4 from a text book, which I have checked several times.



I have a very limited scoop when it comes to numerical programming and I am trying to get to better grips with i, and understand if there more efficient way of coding this but if possible I would like to understand what incorrect within my code, if possible. Any help would be much appreciated.










share|cite|improve this question









$endgroup$








  • 1




    $begingroup$
    What you're asking is not very user friendly. You drop us some code (without any comments... not a good start for programming) and ask us to debug it. At least, it seems that Desmos output is the one of damped oscillator while the original oscillator equation has no damping. You should look in that direction.
    $endgroup$
    – mathcounterexamples.net
    Jan 7 at 20:50










  • $begingroup$
    Your g in k2v = g(...) is only given two arguments. There is a + which should've been a ,. Same with k4v.
    $endgroup$
    – Arthur
    Jan 7 at 20:53












  • $begingroup$
    Also heads up: RK4 is not energy preserving, so even with your error "corrected", you may notice dampening anyway. For this reason, people don't use RK4 for astrophysics simulations.
    $endgroup$
    – Cameron Williams
    Jan 7 at 21:29










  • $begingroup$
    Thank you for the comments, firstly and my appoliges for not putting comments in, that I will take note of.
    $endgroup$
    – james2018
    Jan 7 at 22:15










  • $begingroup$
    @Cameron Williams I am not sure that the OP targets astrophysical calculations... Being myself an ordinary user of Matlab RK4, I can say that, for well behaved ODE,you have to wait a pretty long time before noticing damping...
    $endgroup$
    – Jean Marie
    Jan 11 at 18:25
















0












$begingroup$


So I am having a issue plotting a simply harmonic motion of the form



$$frac{d^2y}{dx^2}+frac{k}{m}y=0$$



Using the RK4 method in matlab. The issue is that my code is not producing the expected plotted and I am not entirely sure, if it my RK4 that is wrong or my actual code that is wrong.



So the solution to the above equation is $y=3cos(frac{k}{m}x)$, if i set k=1 and m=1, I should produce the following graph.



enter image description here



Using Demos



The graph that is being produced by my code is



enter image description here



here is the code which produces the given graph



clc;
clear all;
%Defing intial paramters
%dt=0.001;

h=0.001;
k=1;
m=1;
t=0:h:100;

x=zeros(length(t),1);
v=zeros(length(t),1);
x(1)=3;
v(1)=0;

%defing functions

f=@(t,x,v) v;
g=@(t,x,v) -k/m*x;





for n=1:(length(t)-1);
k1x=f(t(n),x(n),v(n));
k1v=g(t(n),x(n),v(n));

k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k3x=f(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);
k3v=g(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);

x(n+1)=x(n)+h/6*(k1x+2*k2x+2*k3x+k4x);
v(n+1)=v(n)+h/6*(k1v+2*k2v+2*k3v+k4v);


end

plot(t,x);


`
I cannot understand to why my graph is being produced this way, I have gone through the code a couple of time, but cant see to find what is wrong. My only conclusion so far is that it either a variable I have missed or place wrong or my for stamen t is incorrect but cant see that being the case due to, using the equation for RK4 from a text book, which I have checked several times.



I have a very limited scoop when it comes to numerical programming and I am trying to get to better grips with i, and understand if there more efficient way of coding this but if possible I would like to understand what incorrect within my code, if possible. Any help would be much appreciated.










share|cite|improve this question









$endgroup$








  • 1




    $begingroup$
    What you're asking is not very user friendly. You drop us some code (without any comments... not a good start for programming) and ask us to debug it. At least, it seems that Desmos output is the one of damped oscillator while the original oscillator equation has no damping. You should look in that direction.
    $endgroup$
    – mathcounterexamples.net
    Jan 7 at 20:50










  • $begingroup$
    Your g in k2v = g(...) is only given two arguments. There is a + which should've been a ,. Same with k4v.
    $endgroup$
    – Arthur
    Jan 7 at 20:53












  • $begingroup$
    Also heads up: RK4 is not energy preserving, so even with your error "corrected", you may notice dampening anyway. For this reason, people don't use RK4 for astrophysics simulations.
    $endgroup$
    – Cameron Williams
    Jan 7 at 21:29










  • $begingroup$
    Thank you for the comments, firstly and my appoliges for not putting comments in, that I will take note of.
    $endgroup$
    – james2018
    Jan 7 at 22:15










  • $begingroup$
    @Cameron Williams I am not sure that the OP targets astrophysical calculations... Being myself an ordinary user of Matlab RK4, I can say that, for well behaved ODE,you have to wait a pretty long time before noticing damping...
    $endgroup$
    – Jean Marie
    Jan 11 at 18:25














0












0








0





$begingroup$


So I am having a issue plotting a simply harmonic motion of the form



$$frac{d^2y}{dx^2}+frac{k}{m}y=0$$



Using the RK4 method in matlab. The issue is that my code is not producing the expected plotted and I am not entirely sure, if it my RK4 that is wrong or my actual code that is wrong.



So the solution to the above equation is $y=3cos(frac{k}{m}x)$, if i set k=1 and m=1, I should produce the following graph.



enter image description here



Using Demos



The graph that is being produced by my code is



enter image description here



here is the code which produces the given graph



clc;
clear all;
%Defing intial paramters
%dt=0.001;

h=0.001;
k=1;
m=1;
t=0:h:100;

x=zeros(length(t),1);
v=zeros(length(t),1);
x(1)=3;
v(1)=0;

%defing functions

f=@(t,x,v) v;
g=@(t,x,v) -k/m*x;





for n=1:(length(t)-1);
k1x=f(t(n),x(n),v(n));
k1v=g(t(n),x(n),v(n));

k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k3x=f(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);
k3v=g(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);

x(n+1)=x(n)+h/6*(k1x+2*k2x+2*k3x+k4x);
v(n+1)=v(n)+h/6*(k1v+2*k2v+2*k3v+k4v);


end

plot(t,x);


`
I cannot understand to why my graph is being produced this way, I have gone through the code a couple of time, but cant see to find what is wrong. My only conclusion so far is that it either a variable I have missed or place wrong or my for stamen t is incorrect but cant see that being the case due to, using the equation for RK4 from a text book, which I have checked several times.



I have a very limited scoop when it comes to numerical programming and I am trying to get to better grips with i, and understand if there more efficient way of coding this but if possible I would like to understand what incorrect within my code, if possible. Any help would be much appreciated.










share|cite|improve this question









$endgroup$




So I am having a issue plotting a simply harmonic motion of the form



$$frac{d^2y}{dx^2}+frac{k}{m}y=0$$



Using the RK4 method in matlab. The issue is that my code is not producing the expected plotted and I am not entirely sure, if it my RK4 that is wrong or my actual code that is wrong.



So the solution to the above equation is $y=3cos(frac{k}{m}x)$, if i set k=1 and m=1, I should produce the following graph.



enter image description here



Using Demos



The graph that is being produced by my code is



enter image description here



here is the code which produces the given graph



clc;
clear all;
%Defing intial paramters
%dt=0.001;

h=0.001;
k=1;
m=1;
t=0:h:100;

x=zeros(length(t),1);
v=zeros(length(t),1);
x(1)=3;
v(1)=0;

%defing functions

f=@(t,x,v) v;
g=@(t,x,v) -k/m*x;





for n=1:(length(t)-1);
k1x=f(t(n),x(n),v(n));
k1v=g(t(n),x(n),v(n));

k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k3x=f(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);
k3v=g(t(n)+0.5*h,x(n)+0.5*h*k2x,v(n)+0.5*h*k2v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);

x(n+1)=x(n)+h/6*(k1x+2*k2x+2*k3x+k4x);
v(n+1)=v(n)+h/6*(k1v+2*k2v+2*k3v+k4v);


end

plot(t,x);


`
I cannot understand to why my graph is being produced this way, I have gone through the code a couple of time, but cant see to find what is wrong. My only conclusion so far is that it either a variable I have missed or place wrong or my for stamen t is incorrect but cant see that being the case due to, using the equation for RK4 from a text book, which I have checked several times.



I have a very limited scoop when it comes to numerical programming and I am trying to get to better grips with i, and understand if there more efficient way of coding this but if possible I would like to understand what incorrect within my code, if possible. Any help would be much appreciated.







numerical-methods physics matlab






share|cite|improve this question













share|cite|improve this question











share|cite|improve this question




share|cite|improve this question










asked Jan 7 at 20:22









james2018james2018

687




687








  • 1




    $begingroup$
    What you're asking is not very user friendly. You drop us some code (without any comments... not a good start for programming) and ask us to debug it. At least, it seems that Desmos output is the one of damped oscillator while the original oscillator equation has no damping. You should look in that direction.
    $endgroup$
    – mathcounterexamples.net
    Jan 7 at 20:50










  • $begingroup$
    Your g in k2v = g(...) is only given two arguments. There is a + which should've been a ,. Same with k4v.
    $endgroup$
    – Arthur
    Jan 7 at 20:53












  • $begingroup$
    Also heads up: RK4 is not energy preserving, so even with your error "corrected", you may notice dampening anyway. For this reason, people don't use RK4 for astrophysics simulations.
    $endgroup$
    – Cameron Williams
    Jan 7 at 21:29










  • $begingroup$
    Thank you for the comments, firstly and my appoliges for not putting comments in, that I will take note of.
    $endgroup$
    – james2018
    Jan 7 at 22:15










  • $begingroup$
    @Cameron Williams I am not sure that the OP targets astrophysical calculations... Being myself an ordinary user of Matlab RK4, I can say that, for well behaved ODE,you have to wait a pretty long time before noticing damping...
    $endgroup$
    – Jean Marie
    Jan 11 at 18:25














  • 1




    $begingroup$
    What you're asking is not very user friendly. You drop us some code (without any comments... not a good start for programming) and ask us to debug it. At least, it seems that Desmos output is the one of damped oscillator while the original oscillator equation has no damping. You should look in that direction.
    $endgroup$
    – mathcounterexamples.net
    Jan 7 at 20:50










  • $begingroup$
    Your g in k2v = g(...) is only given two arguments. There is a + which should've been a ,. Same with k4v.
    $endgroup$
    – Arthur
    Jan 7 at 20:53












  • $begingroup$
    Also heads up: RK4 is not energy preserving, so even with your error "corrected", you may notice dampening anyway. For this reason, people don't use RK4 for astrophysics simulations.
    $endgroup$
    – Cameron Williams
    Jan 7 at 21:29










  • $begingroup$
    Thank you for the comments, firstly and my appoliges for not putting comments in, that I will take note of.
    $endgroup$
    – james2018
    Jan 7 at 22:15










  • $begingroup$
    @Cameron Williams I am not sure that the OP targets astrophysical calculations... Being myself an ordinary user of Matlab RK4, I can say that, for well behaved ODE,you have to wait a pretty long time before noticing damping...
    $endgroup$
    – Jean Marie
    Jan 11 at 18:25








1




1




$begingroup$
What you're asking is not very user friendly. You drop us some code (without any comments... not a good start for programming) and ask us to debug it. At least, it seems that Desmos output is the one of damped oscillator while the original oscillator equation has no damping. You should look in that direction.
$endgroup$
– mathcounterexamples.net
Jan 7 at 20:50




$begingroup$
What you're asking is not very user friendly. You drop us some code (without any comments... not a good start for programming) and ask us to debug it. At least, it seems that Desmos output is the one of damped oscillator while the original oscillator equation has no damping. You should look in that direction.
$endgroup$
– mathcounterexamples.net
Jan 7 at 20:50












$begingroup$
Your g in k2v = g(...) is only given two arguments. There is a + which should've been a ,. Same with k4v.
$endgroup$
– Arthur
Jan 7 at 20:53






$begingroup$
Your g in k2v = g(...) is only given two arguments. There is a + which should've been a ,. Same with k4v.
$endgroup$
– Arthur
Jan 7 at 20:53














$begingroup$
Also heads up: RK4 is not energy preserving, so even with your error "corrected", you may notice dampening anyway. For this reason, people don't use RK4 for astrophysics simulations.
$endgroup$
– Cameron Williams
Jan 7 at 21:29




$begingroup$
Also heads up: RK4 is not energy preserving, so even with your error "corrected", you may notice dampening anyway. For this reason, people don't use RK4 for astrophysics simulations.
$endgroup$
– Cameron Williams
Jan 7 at 21:29












$begingroup$
Thank you for the comments, firstly and my appoliges for not putting comments in, that I will take note of.
$endgroup$
– james2018
Jan 7 at 22:15




$begingroup$
Thank you for the comments, firstly and my appoliges for not putting comments in, that I will take note of.
$endgroup$
– james2018
Jan 7 at 22:15












$begingroup$
@Cameron Williams I am not sure that the OP targets astrophysical calculations... Being myself an ordinary user of Matlab RK4, I can say that, for well behaved ODE,you have to wait a pretty long time before noticing damping...
$endgroup$
– Jean Marie
Jan 11 at 18:25




$begingroup$
@Cameron Williams I am not sure that the OP targets astrophysical calculations... Being myself an ordinary user of Matlab RK4, I can say that, for well behaved ODE,you have to wait a pretty long time before noticing damping...
$endgroup$
– Jean Marie
Jan 11 at 18:25










1 Answer
1






active

oldest

votes


















2












$begingroup$

I do not know why that does not produce an error of insufficient arguments. You have two identical typos where instead of a comma there is a plus in the computation of k2v and k4v.



k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);





share|cite|improve this answer









$endgroup$













    Your Answer





    StackExchange.ifUsing("editor", function () {
    return StackExchange.using("mathjaxEditing", function () {
    StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
    StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
    });
    });
    }, "mathjax-editing");

    StackExchange.ready(function() {
    var channelOptions = {
    tags: "".split(" "),
    id: "69"
    };
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function() {
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled) {
    StackExchange.using("snippets", function() {
    createEditor();
    });
    }
    else {
    createEditor();
    }
    });

    function createEditor() {
    StackExchange.prepareEditor({
    heartbeatType: 'answer',
    autoActivateHeartbeat: false,
    convertImagesToLinks: true,
    noModals: true,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: 10,
    bindNavPrevention: true,
    postfix: "",
    imageUploader: {
    brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
    contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
    allowUrls: true
    },
    noCode: true, onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    });


    }
    });














    draft saved

    draft discarded


















    StackExchange.ready(
    function () {
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3065455%2fsimple-harmonic-motion-rk4-method-matlab-plot%23new-answer', 'question_page');
    }
    );

    Post as a guest















    Required, but never shown

























    1 Answer
    1






    active

    oldest

    votes








    1 Answer
    1






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes









    2












    $begingroup$

    I do not know why that does not produce an error of insufficient arguments. You have two identical typos where instead of a comma there is a plus in the computation of k2v and k4v.



    k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
    k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

    k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
    k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);





    share|cite|improve this answer









    $endgroup$


















      2












      $begingroup$

      I do not know why that does not produce an error of insufficient arguments. You have two identical typos where instead of a comma there is a plus in the computation of k2v and k4v.



      k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
      k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

      k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
      k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);





      share|cite|improve this answer









      $endgroup$
















        2












        2








        2





        $begingroup$

        I do not know why that does not produce an error of insufficient arguments. You have two identical typos where instead of a comma there is a plus in the computation of k2v and k4v.



        k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
        k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

        k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
        k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);





        share|cite|improve this answer









        $endgroup$



        I do not know why that does not produce an error of insufficient arguments. You have two identical typos where instead of a comma there is a plus in the computation of k2v and k4v.



        k2x=f(t(n)+0.5*h,x(n)+0.5*h*k1x,v(n)+0.5*h*k1v);
        k2v=g(t(n)+0.5*h,x(n)+0.5*h*k1x+v(n)+0.5*h*k1v);

        k4x=f(t(n)+h,x(n)+h*k3x,v(n)+h*k3v);
        k4v=g(t(n)+h,x(n)+h*k3x+v(n)+h*k3v);






        share|cite|improve this answer












        share|cite|improve this answer



        share|cite|improve this answer










        answered Jan 7 at 21:25









        LutzLLutzL

        57.3k42054




        57.3k42054






























            draft saved

            draft discarded




















































            Thanks for contributing an answer to Mathematics Stack Exchange!


            • Please be sure to answer the question. Provide details and share your research!

            But avoid



            • Asking for help, clarification, or responding to other answers.

            • Making statements based on opinion; back them up with references or personal experience.


            Use MathJax to format equations. MathJax reference.


            To learn more, see our tips on writing great answers.




            draft saved


            draft discarded














            StackExchange.ready(
            function () {
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3065455%2fsimple-harmonic-motion-rk4-method-matlab-plot%23new-answer', 'question_page');
            }
            );

            Post as a guest















            Required, but never shown





















































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown

































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown







            Popular posts from this blog

            1300-talet

            1300-talet

            Display a custom attribute below product name in the front-end Magento 1.9.3.8