How to plot a nonlinear equation using FSOLVE in MATLAB

Hello

I have a nonlinear equation like this:

in which:

I want to plot vs. k

obviously this is a nonlinear equation. How should I use FSOLVE to have this plot?

Answer 1

[John D'Errico]

You have an implicit function. So a relationship involving two variables, where there is no way to resolve one as a function of the other directly. You cannot write ky = f(k) only.

We might use the symbolic toolbox. But there is no direct need to do so. Just because you don't know the value of two variables, does not mean you must have them be symbolic unknowns. A simple function will suffice. However, you CANNOT easily use fsolve. Why not? Because fsolve will give you only ONE solution!

kp = 4500;
h = 1e-3;
d = 7.5e-3;
fun = @(k,ky) ky./k.*tan(ky*h)+(1-(k.^2-ky.^2)./(kp^2+k.^2-ky.^2)).*tan(k*d)+(k.^2-ky.^2)./(kp^2+k.^2-ky.^2).*sqrt(kp^2-k.^2)./k.*tan(sqrt(kp^2-ky.^2)*d);

See that fun has been writtten to allow for vectorized computation. So we can pass in vectors or arrays of numbers for k and ky. To that end, I used the .*, ./, and .^ operators. You can use the non-dotted operators where a scalar operation is involved. So ky*h is sufficient, and fully vectorized.

Now just plot it using fimplicit. I had to play around with the axis limits on ky, to see the full behavior of this relationship.

fimplicit(fun,[0 628 0 5000])
xlabel k
ylabel ky

fimplicit shows us there are MANY solutions for any given value of k. In fact, it looks at first like there may be infinitely many of them. But then we see the plot seems to cut off at the top end, around ky at 4500. At this point, you might look more carefully at what happens at ky==4500, and it becomes clear that is indeed the upper limit. We would probably begin to have complex solutions only above that point for ky.

The problem with using fsolve on this, is that fsolve will not tell you the story. And fsolve will give you just one solution, that WILL depend on the starting value. Therefore, fsolve will tell you an incomplete story, and very probably a confusing one.

For example, suppose you pick some value of k. Say k==300. Substiture k == 300. Now, you could use fsolve (or zfero) to then solve for ONE value of ky. But you could solve for roughly 20 values of ky at any given k.

Anyway, having now resolved what this implicit function looks like, we could now refine the limits. Or, we could spend some time in investigating just why it is the function looks as it does. But you need to fully visualize what happens before that is at all possible. For example, now, try it out. For kicks, I'll use the symbolic TB here.

syms k ky
sfun = ky./k.*tan(ky*h)+(1-(k.^2-ky.^2)./(kp^2+k.^2-ky.^2)).*tan(k*d)+(k.^2-ky.^2)./(kp^2+k.^2-ky.^2).*sqrt(kp^2-k.^2)./k.*tan(sqrt(kp^2-ky.^2)*d);
pretty(subs(fun,k,100))
                                                                     /                     2  \
      /  ky  \                                                       | 3 sqrt(20250000 - ky ) |    2
ky tan| ---- |            /     2              \   sqrt(20240000) tan| ---------------------- | (ky  - 10000)
      \ 1000 /      / 3 \ |   ky  - 10000      |                     \           400          /
-------------- - tan| - | | -------------- - 1 | + ----------------------------------------------------------
      100           \ 4 / |   2                |                         2
                          \ ky  - 20260000     /                      (ky  - 20260000) 100

UGH. It seems clear there will be no direct solution. Can solve do some magic anyway?

solve(subs(fun,k,100))
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve. 
> In solve (line 304) 
ans =
1031.1951300726862936928503667825

So solve did indeed give up. When it does so, it passes things off to vpasolve. Again though, this is a numerical solver. So it will find only one solution. If we give vpasolve a different start point, it will tell a different story, just as fzero or fsolve would have done.

vpasolve(subs(fun,k,100),2000)
ans =
2029.1205275928410282293585811077

All using the symbolic toolbox did for me was to produce more accurate solutions. Lets return to fun. I'll pick just one value of k, say at k == 100. Now look at the various solutions. fplot will be sufficient here.

fplot(@(ky) fun(100,ky),[0 4500])
yline(0);
xlabel ky
title 'fun, at a fixed value of k==100'
ylabel 'fun(100,ky)'

Better now. We can see that the horizontal lines in the fimplicit plot were not true solutions, but where it saw the function cross between singularities. So those horizontal lines in the fimplicit plots were spurious solutions. fplot shows them as dotted vertical lines here. The horizontal reference line indicates where the function crosses zero. So, at k == 100, we should expect to see the first solution at ky just a bit above 1000, and then a second solution just over 2000, etc. There is one other point of interest perhaps, and that might be to investigate what happens around k==210 or so. How would we do that? Simplest seems to look at the symbolic version of this implicit function.

pretty(sfun)
                                                             /                     2  \
      /  ky  \                                               | 3 sqrt(20250000 - ky ) |                  2    2     2
ky tan| ---- |              /        2     2          \   tan| ---------------------- | sqrt(20250000 - k ) (k  - ky )
      \ 1000 /      / 3 k \ |       k  - ky           |      \           400          /
-------------- - tan| --- | | ------------------- - 1 | + ------------------------------------------------------------
       k            \ 400 / |  2     2                |                          2     2
                            \ k  - ky  + 20250000     /                      k (k  - ky  + 20250000)

Anything strike you? When I look at that mess, I am looking for something that involves strictly k. I hope what stands out is the sub-term tan(3*k/400). (Yes, it seems to be a real mess here, because the site is not doing a good job of formatting this specific content. In the MATLAB command window, things were a bit more clear.)

The tangent function goes through a singularity for tan(theta) when theta = pi/2, and again at theta = 3*pi/2. Essntially at pi*(n + 1/2), for integer values of n. This should tell us where the singularity should arise for k. Just solve for k, such that

3*k/400 == pi*(n + 1/2)

That gives us a singularity where things get nasty, at:

k = 400/3*pi*(n + 1/2)

The first such singularity happens at

400/3*pi/2
ans =
       209.44

And that seems to be pretty much the complete story.

Answer 2

[Torsten]

K = 1:627;
kp = 4500;
h = 1.0e-3;
d = 7.5e-3;
ky0 = 1.0;
Ky = zeros(numel(K),1);
%fKy = zeros(numel(K),1);
for i = 1:numel(K)
  k = K(i);
  fun = @(ky) ky/k*tan(ky*h)+(1-(k^2-ky^2)/(kp^2+k^2-ky^2))*tan(k*d)+(k^2-ky^2)/(kp^2+k^2-ky^2)*sqrt(kp^2-k^2)/k*tan(sqrt(kp^2-ky^2)*d);
  Ky(i) = fzero(fun,ky0);
  %fKy(i) = fun(Ky(i));
  ky0 = Ky(i);
end
plot(K,Ky) 
%plot(K,fKy)

Answer 3

[mohammaddmt]

Thank you both so much for your answers.

let me put the problem this way:

the plot should be something like this:

in which ky is plotted vs. f.

I asked my prof. and he said you should consider both real and imaginary parts of ky.

so imaginary part of fun should be like blue line and the real part should be like red line.

how could I possibly break it down to real and imaginary parts?

[John D'Errico]

Please stop adding answers just to make a comment.

[John D'Errico]

Anyway, I just spent what, well over an hour writing a solution to the problem that you posed, i.e., how to use fsolve to plot the solutions of an equation. Since fsolve only returns real solutions, that question was actually irrelevant, as was much of my solution as written.

Is there a good reason why you did not tell us what you really wanted to learn, and thus save me a great deal of time in answering the wrong question?

[mohammaddmt]

I'm really really sorry I didn't mention this, since I just myself asked my Prof. to help me in this regard and he showed me this figure. but anyway, I'm deeply thankful to you. I learned so much from your asnwer.

I should probably figure that out on my own based on your precious answer.

Thank you so much for the effort you made.

If only you could help me with the hint that my Prof. told, I would be very thankful to you.

should I ask a new question?

Again I'm sorry for what occured.

BR

Answer 4

[Mahmoud Abdelaziz]

Hi there,

I see an excellent and helpfull answer.

I have a question. Can you help me.

I want to plot some nonlinear inequalities as follows:

R1>0,R2<0,R3>0,R4>0,R5>0,R6>0, where

R1= x*y+2, R2=x-1, R3=x^2*y-x*y-1, R4=x^2*y-x*y+2, R5=x*y+3, R6=x^3*y-4*x^2*y+3*x^2+4*x*y-8*x+8.

If not possible to plot these inequalities and find the region of solutions, Can plot them as equations?

Thank youfor your help.

M. A. M. Abdelaziz.

[John D'Errico]

Ask this as a separate question. What you have posted is not an answer, but a question.

By the way, what you have shown are NOT inequalities, but equalities. There is a difference, a big one. There is no region of solution to what you have posted, where they are equality statements.

Regardless, if you want help, post this as a QUESTION.