How the increment size would affect solving an equation ?

I have a function, and this function should calculate simple equation in a range from 0 to pi/2; the main issue is that I would like to have accurate results, and this require small increment. Here is the function:

function [Sg_wk]=MM(M1,delta)
Y=1.4;
Sg=0:0.001:pi/2;
Del=zeros(1,length(Sg));
    for i=1:1:length(Sg)
    Del(i)=atan((2*cot(Sg(i))*(M1^2*sin(Sg(i))^2-1))/( 2+M1^2*(Y+cos(2*Sg(i))))); 
    end
Sg_wk=interp1(rad2deg(Del), rad2deg(Sg), delta, 'linear');
end

If I run this code with >> MM(2,4) ans = 48.6861

The correct answer is 33.39024 If I use the same code with smaller increment (Sg=0:0.1:pi/2), it provides: >> MM(2,4)

ans = 33.4677

much better, but does not that contradict the logic ???

Also, Del now will be calculated, how can extract only the positive values out of the whole matrix ?

Answer 1

[John D'Errico]

Lets make this easy.

M1 = 2;
delta = 4;
Y=1.4;
D = @(Sg) atan((2*cot(Sg).*(M1^2*sin(Sg).^2-1))./( 2+M1^2*(Y+cos(2*Sg))));

Sg=0:0.001:pi/2;
Del = D(Sg);

Now, plot what you are trying to use for interpolation. The INDEPENDENT variable is essentially Del.

plot(rad2deg(Del), rad2deg(Sg))
grid on

So, if you try to interpolate that curve from that data for delta between 0 and 20, there are TWO values. Your function is multi-valued.

You CANNOT use interp1 to interpolate that data. Interp1 assumes a single valued function, something you do not have.

Can you solve the problem? Using syms...

syms Sg
M1 = 2;
Y=1.4;

vpasolve(atan((2*cot(Sg).*(M1^2*sin(Sg).^2-1))./( 2+M1^2*(Y+cos(2*Sg)))) == deg2rad(4),Sg,.3)
ans =
0.58276955173511339987159569253578

rad2deg(double(ans))
ans =
      33.39

Now, try starting it out so it will find the other solution.

vpasolve(atan((2*cot(Sg).*(M1^2*sin(Sg).^2-1))./( 2+M1^2*(Y+cos(2*Sg)))) == deg2rad(4),Sg,1.5)
ans =
1.5285992045896074581678383895601

rad2deg(double(ans))
ans =
       87.582

So, there are two solutions, one around 34 degrees, the other around 87 degrees.

I could have done this using fzero too.

[Saad Moalwi]

Yes, you are right, for some deltas (x-axis) there will be two sigmas (y-axis). So what command instead of interp1 will help to find the two sigmas. (My code is lengthy, and simple change will yield to rebuild the whole thing, thus I can only change the interp1 command).

[John D'Errico]

You cannot use interpolation directly in this way. Interpolation predicts a SINGLE valued function. It makes no sense at all to use it for a multi-valued problem.

You have two choices that I see.

1. Use Doug Schwarz's intersections code. You will pass it the curve that you have as one curve. Then pass it a curve defined by TWO points. Intersections will find the points of intersection, if any.

intersections(rad2deg(Sg), rad2deg(Del),[0 90],[delta,delta])
ans =
      33.39
     87.582

2. Alternatively, you can download my SLM toolbox , also from the File Exchange. This would give a slightly more accurate solution, based on use of a spline interpolant, whereas intersections will use linear interpolation.

pp = pchip(rad2deg(Sg), rad2deg(Del));
slmsolve(pp,4)
ans =
        33.39       87.582

Either will be entirely adequate here. As you can see, they are identical to within 5 significant digits on this problem.