Stuck with matlab

[dugironbeard]

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.

[uart]

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.

[dugironbeard]

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?

[uart]

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?

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.

[dugironbeard]

Alright,

Well im going back in tomorrow. Ill google a bit how to make a script for this shooting method.

[dugironbeard]

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);