Usage of fzero together with ode45 solution

I have a differential equation which I solve with ode45 like this:

[x,u]=ode45(@(x,u) diffljud2(x,u,v), [0,151900], [2000, tand(v)]); fun=u(:,1); %z(x)

The task is to find values for -10<v<14 that makes z(x) close to 2500 when x=151900 I tried like this:

func=@(v) u(:,1)-2500; v0=-10; v=fzero(func, v0)

I get the error "undefined function or variable v", though I'm not sure that I'm creating the function in the right way..

[Star Strider]

In order to attempt to help, we need as much of your code as necessary to integrate your ‘diffljud2’ ODE function, including ‘u’ and ‘v’.

[Tomas Hermansson]

Okay, here's the function: function diffekv=diffljud2(x,u,v)

p1=-22.0335979129617; p2=17.6231123495468; p3=274.5871624167000; p4=0.7457214382711;

cz=(4800+p1+p2*(u(1)/1000)+p3*exp(-p4*(u(1)/1000)))^3; %

czprim=p2/1000-p3*(p4/1000)*exp(-p4*(u(1)/1000));

q0=(funkljudvagor2(2000)/cosd(v))^2; %funkljudvagor(2000) is a constant

diffekv=[u(2); -q0*(czprim/cz)];

end

[Jan]

@Tomas: Please use the "{} Code" button to format your code. Currently it is hard to read.

If you post the complete error message by Copy&Paste, would could see which line causes the problem.

[Star Strider]

I believe we still need ‘v’.

I see that ‘v’ is defined in the fzero result, however it is a function of ‘u’ that is the result of integrating your differential equation.

Do you need to do the fzero call inside ‘diffljud2’?

Please go into some detail as to what you want to do. I do not understand your code.

Answer 1

[Jan]

Is "v" defined before the line:

[x,u] = ode45(@(x,u) diffljud2(x,u,v), [0,151900], [2000, tand(v)])

?

What is the purpose of

 func = @(v) u(:,1)-2500;
 v0   = -10;
 v    = fzero(func, v0)

The anonymous function "func" takes the variable "v" as input, but uses the value of "u"? Do you want to find the row of "u", which is nearest to 2500? Then:

[value, index] = min(abs(u(:,1) - 2500))

Answer 2

[Tomas Hermansson]

Okay, probably I should post my previous code to clarify what I'm really looking for. This is the task

Now suppose that a sound source at a depth of 2000 feet transmits to a receiver 25 miles away at a depth of 2500 feet. The above calculation shows that one of the rays from the source to the receiver leaves the source at an angle close to 7.9 degrees. Because of the nonlinearity of the equation there may be other rays leaving at different angles that reach the same receiver. Run your program for β0 in the range from −10 up to 14 degrees, plot the ray paths and print a table of the values z(xf) (xf=151900 feet=25 nautical miles). We are interested in finding values of β0 for which z(xf ) = 2500. Use an efficient algorithm to determine the rays which pass through the receiver. Discuss the accuracy of your results.

I've plotted the rays from -10 to 14 with this code:

for v=-10:0.1:14

    tolerans=odeset('RelTol',1e-6);

    [x,u]=ode45(@(x,u)diffljud2(x,u,v), [0,151900], [2000, tand(v)],tolerans);

    plot(x,u(:,1),'-')
    hold on
    grid on
    axis ij
end

That gives me a nice plot of the rays. From this you can see that u(:,1)=z(x). I think my biggest trouble is how to use u(:,1) with fzero to find the angles that the assignment asks for.

Thanks for trying to help!

[Star Strider]

I just encountered another problem:

Undefined function or variable 'funkljudvagor2'.

      q0=(funkljudvagor2(2000)/cosd(v))^2; %funkljudvagor(2000) is a constant

[Tomas Hermansson]

Oh, just change funkljudvagor(2000) to 4.875007837540556e+03 instead and it should work

[Jan]

@Tomas: Please do not post details of the question in the section for answer. This is confusing.

"u" is a vector. It is simply not sufficient as input for fzero, which requires a function. Vectors are not functions and vice versa.

Note that this question has been discussed repeatedly in this and other Matlab forums. So a web search might help you.

If you want to use fzero, the function must include the integration and the optimal v is determined as root. Lookfor "single shooting" methods for more details.

Answer 3

[Star Strider]

The interp1 function may be more useful than fzero here. The output vector of the result I believe you want is ‘vx’.

See if this does what you want:

function diffekv=diffljud2(x,u,v)
p1=-22.0335979129617;
p2=17.6231123495468;
p3=274.5871624167000;
p4=0.7457214382711;
cz=(4800+p1+p2*(u(1)/1000)+p3*exp(-p4*(u(1)/1000)))^3; %
czprim=p2/1000-p3*(p4/1000)*exp(-p4*(u(1)/1000));
q0=(4.875007837540556e+03/cosd(v))^2; %funkljudvagor(2000) is a constant
diffekv=[u(2); -q0*(czprim/cz)];
end

v=-10:0.1:14;
for k1=1:length(v)

      tolerans=odeset('RelTol',1e-6);

      [x,u]=ode45(@(x,u)diffljud2(x,u,v(k1)), [0,151900], [2000, tand(v(k1))],tolerans);

      vx(k1) = interp1(u(:,1), x, 2500, 'linear','extrap');

      plot(x,u(:,1),'-')
      hold on
      grid on
      axis ij
  end

figure(2)
plot(v, vx)
grid

You will likely need to experiment, since I still have no idea what you are doing.

This code runs. You must determine if it gives the result you want.