Simple harmonic motion, RK4 method, matlab plot
$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.
Using Demos
The graph that is being produced by my code is
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
asked Jan 7 at 20:22
james2018james2018
827
$endgroup$
add a comment |
$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.
Using Demos
The graph that is being produced by my code is
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
asked Jan 7 at 20:22
james2018james2018
827
$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$
Yourgink2v = g(...)is only given two arguments. There is a+which should've been a,. Same withk4v.
$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
add a comment |
$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.
Using Demos
The graph that is being produced by my code is
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
asked Jan 7 at 20:22
james2018james2018
827
$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.
Using Demos
The graph that is being produced by my code is
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
numerical-methods physics matlab
asked Jan 7 at 20:22
james2018james2018
827
asked Jan 7 at 20:22
james2018james2018
827
asked Jan 7 at 20:22
james2018james2018
827
asked Jan 7 at 20:22
james2018james2018
827
asked Jan 7 at 20:22
james2018james2018
827
827
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$
Yourgink2v = g(...)is only given two arguments. There is a+which should've been a,. Same withk4v.
$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
add a comment |
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$
Yourgink2v = g(...)is only given two arguments. There is a+which should've been a,. Same withk4v.
$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
add a comment |
1 Answer
1
active
oldest
votes
$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);
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
$endgroup$
add a comment |
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
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
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
$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);
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
$endgroup$
add a comment |
$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);
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
$endgroup$
add a comment |
$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);
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
$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);
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
answered Jan 7 at 21:25
LutzLLutzL
58.7k42055
58.7k42055
add a comment |
add a comment |
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.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
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
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
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
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
gink2v = g(...)is only given two arguments. There is a+which should've been a,. Same withk4v.$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