1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Stuck with matlab

  1. Oct 27, 2009 #1
    1. The problem statement, all variables and given/known data

    dx(1)/dt=-0.1x(1)-x(2)
    dx(2)/dt=x(1)-0.1x(2)

    x1|t=0 = 1 , x2|t=0 =0

    First we had to define the magnitude (q) of the velocity vector at t = 10.

    The follow up question is:

    The horizontal component x1|t=0 of the starting position is not known very accurately. However at t = 10 an accurate measurement of q=0.45 is known.

    Determine the vvalue of x1|t=0 that leads to q=0.45

    2. Relevant equations



    3. The attempt at a solution


    To find the initial q I made the following files in matlab.
    clear,clc

    x0=[1 0];
    tspan=[0 1.01 10];
    [t,x]=ode45('partb',tspan,x0);

    x =

    1.0000 0
    0.4807 0.7656
    -0.3083 -0.2002


    function xp = partb(t,x)
    t=10
    x(1)=-0.3083
    x(2)=-0.2002
    xp(1)=(-0.1*x(1)-x(2));
    xp(2)=(x(1)-0.1*x(2));
    disp (xp(1))
    disp (xp(2))
    q=((xp(1)^2)+(xp(2)^2))^0.5
    end

    q = 0.3694

    Now the problem is I dont know how to work backwards, I know I have to do something with fzero but I have to many variables. Any help would be greatly appreciated.
     
  2. jcsd
  3. Oct 27, 2009 #2

    uart

    User Avatar
    Science Advisor

    It's called the "shooting method". You need to perturbate the initial conditions (IC's) to search for the correct IC to get x1(10)=0.45.

    Start by finding one IC that gives x1(10)<0.45 and one other that gives x1(10)>0.45 then just narrow it down using a "root search" algorithm of your choice, though for simplicity I recommend using "halving the interval".

    Note that you have already found one IC that leads to x1(10)<0.45. So now just perturbate x1(0) to find an IC that leads to x1(10)>0.45 and you're good to start halving the interval.
     
  4. Oct 27, 2009 #3
    Thanks a lot for your reply uart.

    Alright so first Im gonna make a IC greater than 0.45. And than I have to write a root search algorithm.

    Secondly really stupid maybe (English is not my native language), what is perturbate?
     
  5. Oct 27, 2009 #4

    uart

    User Avatar
    Science Advisor

    It was probably not quite the right word. I wish I could claim that English was not my native language. :redface: A perturbation in this context is just a "change in" or a disturbance in the initial condition. In other words just try x1(0) either a bit more than one or a bit less than one and see which way x1(10) moves. Once you get it moving in the correct direction then keep going until you get to x1(10) > 0.45, then you are ready to start halving the interval.

    If you like you can initially just do the "halving the interval" manually, making a run of the DE simulation, calculating the new IC (halfway between the current two closest that are located on alternate sides of x1(10)=0.45) and then taking another run etc. You should be able to get an approximate solution fairly quickly like this. One you get the hang of what's happening then you can write some code to implement this "shooting method" automatically.
     
  6. Oct 27, 2009 #5
    Alright,

    Well im going back in tomorrow. Ill google a bit how to make a script for this shooting method.
     
  7. Oct 28, 2009 #6
    Right, next day.

    I have found the initial values of between 1.15 and 1.25. Now to be honest I understand the intention of the shooting method but I have no clue whatsoever how to write it into a matlab code. Problem is I have two different files, one with the initial ode45 and the next one that takes it and calculates a velocity vector.

    What I have now and not working:
    clear,clc
    a=1.20:0.0001:1.25
    x0=[a 0];
    tspan=[0 1.01 10];
    [t,x]=ode45('partb',tspan,x0);
     
    Last edited: Oct 28, 2009
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Stuck with matlab
  1. Stuck with mechanics (Replies: 0)

  2. Grades.java stuck (Replies: 1)

Loading...