Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Diff EQ help

  1. Jul 6, 2005 #1
    I am new to matlab and need help with a problem. we have tried
    several different methods and each time cant make sense of what is going on wrong

    " A mass spring motion is governed by the ordinary differential
    equation m(dx^2/dt^2) + b(dx/dt) + kx=F(t) , where m is the mass, b
    is the damping constant, k is the spring constant, and F(t) is the
    external force."
    Assume the following m=1 kg, k= 16 N/m, and F(t)=0, x(0)=0, x'(0)=0.
    Solve the ODE on the interval [0,20] for the following cases.

    a). b=0 -no damping
    b). b=2- under damping
    c). b=8 - critical damping
    d). b=10 - over damping

    Set the damping constant equal to 0, b=0, and consider an agiing spring whose spring constant depends on time as follows


    Predict the motion of the mass, i.e. show time plots and phase plots and discuss your results for the following cases.

    a. L=0
    b. L=0.1
    c. L=0.3

    Set b = 1/5 k=1/5 and assume a forcing F(t) = coswt,

    a. Use ODE45-solver to obtain the solution curves for several values of w between .5 and 2. Plot the solutions and estimate the amplitude A of the steady response in each case.

    b. Using the data from part (a) plot the graph of A vs. w. For what frequency W is the amplitude greatest?
  2. jcsd
  3. Jul 6, 2005 #2


    User Avatar
    Science Advisor
    Homework Helper

    Hey Gob, you got to push that spring dude. Think about it: You have a mass on a spring sitting there. It has no external force, you don't displace it, and you don't give it an initial velocity. It's not going anywhere, i.e. the only solution is y(x)=0. Are you sure you have the initial conditions right? How about just a little push, say y'(0)=0.1. That'll get it going. :smile:
  4. Jul 6, 2005 #3
    your right, the intial conditions are x(0)=1, x'(0)=0

    the first answer is solved.....we are having problems with #2

    or m file looks something like this

    function out=spring2(t,u)
    global l m b k
    out(1)= u(2);
    out(2)= 16 exp(-l*t)

    trying to run it it gives us u is undefined and we cant figure out what is wrong, and we cant even get to the global part to where we can sub in our different l values
  5. Jul 6, 2005 #4


    User Avatar
    Science Advisor
    Homework Helper

    Sorry, can' help with MatLab. I use Mathematica.
  6. Jul 6, 2005 #5


    User Avatar
    Science Advisor
    Homework Helper

    You know Gob, it makes sense: as you increase the value of L, the rate at which the spring constant decreases, increases right for:


    and with no dampening, as the spring "constant" decreases, the oscillations will get larger, quicker for increasing L. That's what it looks like to me anyway. Can't you just do a numerical solve of the three ODEs with the three values of L in MatLab and plot the results?. In Mathematica it's a breeze.
  7. Jul 6, 2005 #6
    matlab blows a big one, heck i even like maple more....this thing is killin me and we are about to go nuts
  8. Jul 6, 2005 #7

    I don't really understand the matlab code you provided. Is that your function for the ode solver? If you could post anything else you have or more specific problems it may help.

    A quick starter for you
    Code (Text):

    % PF for gob71
    % x(t) is position, y(t)=dx/dt=velocity
    clear all;

    dt = 0.0005; % Time step size
    t = 0:dt:20; % Vector of time values

    % Constants & Parameters
    m = 1; k = 16; b = 0; w = 0; F = 0; L = 0.0;
    y(1)=0; yold = y(1); % Initial Velocity
    x(1)=1; xold = x(1); % Initial Position

    % Real men write their own ODE solvers
    for i = 2:length(t)
        k = exp(-L*t(i));
       % F = cos(w*i);
        x(i) = (m*xold/dt^2 + b*xold/dt + m*yold/dt + F) / (m/dt^2 + b/dt + k);
        y(i) = (x(i)-x(i-1))/dt;
        xold = x(i); yold = y(i);

    A few quick checks for your solution.
    x'' + b/m*x' + k/m*x = F/m
    If b=F=0 and m=1, then when k=16, x(t)=cos(4t), and when k=exp(0*t), x(t)=cos(1*t)
  9. Jul 7, 2005 #8


    User Avatar
    Science Advisor
    Homework Helper

    Why, when you can do this: 1,2,3, done . . .
    Code (Text):

    sol1 = NDSolve[{x''[t] + 16 e^{-L t} x[t] == 0,
          x[0] == 1, x'[0] == 0}, x, {t, 0, 20}]
    f[t_] := Evaluate[x[t] /. Flatten[sol1]];
    Plot[f[u], {u, 0, 20}]
    The plots are for L=0.1, 0.2, 0.3 respectively.

    Attached Files:

    • 01.JPG
      File size:
      8.5 KB
    • 02.JPG
      File size:
      7.4 KB
    • 03.JPG
      File size:
      5.9 KB
    Last edited: Jul 7, 2005
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook