Equation solution is not found for all array element

Hi I have an equation that includes just one variable. With vpasolve command, it is just solved until 105th array value. But array length is 298x1. what is the wrong? and is there any suggestion? (my opinion is, when i equal to 106 and higher, program can't find real solution, but it should be)

clear all, clc, format shortG, close all
name="XXX";
S=load(name+".m");
F=S(:,1);
Vd=S(:,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms t_bulk Te n hw me t_p eps_inf eps_s w
% consts
e=1.602E-19;
h=6.626E-34;
h_p=1.05E-34;
eps_0=8.854E-12;
kb=1.38e-23;
T_lattice=300;
m0=9.1E-31;
thickness=5E-5;
% Parameters
n=(2.59E14)/thickness;
hw=0.034*e;
me=0.057*m0;
t_p=160E-15;
eps_inf=11.46*eps_0;
eps_s=13.94*eps_0;
w=(hw/h_p);
Nc=(2*((2*pi*me*kb*T_lattice)/(h^2))^1.5)*1E-6;
K0=(hw/(2*kb*Te));
t0=1/(((e^2)*(w)/(2*pi*h_p))*sqrt(me/(2*hw))*((1/(eps_inf))-(1/(eps_s))));
t_bulk=(t0+t_p*((n*kb*Te)/(2*Nc*hw))*(1-exp(-hw/(kb*Te))));
n0=(1/(exp(hw/(kb*T_lattice))-1));
Pexp=(e.*F*1000.*Vd);
Ptheory1=((sqrt(2/pi)*(hw/t_bulk)*((n0+1)*exp(-hw/(kb*Te))-n0)*(sqrt(hw/(2*kb*Te)))*exp(hw/(2*kb*Te))*besseli(0,K0)));
Ptheory1=simplify(Ptheory1);

for i=1:length(Pexp)
    T1(i,1)=double(vpasolve(Ptheory1==Pexp(i,1),Te, 1));
end

And error is below: Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 0-by-1.

Answer 1

[John D'Errico]

For SOME of these sub-problems, no solution was found.

I would suggest that when the result from vpasolve is empty, you could assign a NaN for those cases.

for i=1:length(Pexp)
  result = double(vpasolve(Ptheory1==Pexp(i,1),Te, 1));
  if isempty(result)
    T1(i,1) = NaN;
  else
    T1(i,1) = result;
  end
end

[Mustafa]

Yes, after 105th cell, turned as NaN like you wrote into the code. But should be solution. Also there is no big difference between 105th cell and 106th Pexp cell (differency is around 1E-2).

Is there any other way like iteration solution. Ptheory1 is an equation and Pexp includes values. How do they match with appropriate Te values?

[Torsten]

If the values of Pexp are similar from call to call, take the solution of the last call as initial guess of the next call:

T1ini = 1;
for i=1:length(Pexp)
    result = double(vpasolve(Ptheory1==Pexp(i,1),Te, T1ini));
    if isempty(result)
        T1(i,1) = NaN;
    else
        T1(i,1) = result;
        T1ini = result;
    end
end

[Mustafa]

Hi Torsten,

That was really nice approach. thank for it.

But solution still turns NAN after 105th cell. I enter some Te values higher than 105th cell the program found, Ptheory function works.

I think I should accept values in range of error limits, not exact solution (I don't know why it did not work).

How do I write the code so as to accept in error limit? also I guess the error limit will change after in every found value by program because Pexp value change.

Could someone help me on this?

[John D'Errico]

Good for theory. There MAY be a solution. There may be no solution that vpasolve could find. Or you may be wrong about the theory.

For example, find the real solution to

exp(x) - a == 0

As a changes, by only a tiny amount, from eps to -eps, suddenly, no real solution will be found.

But you cannot tell vpasolve to return a solution within some error limit. vpasolve is not a tool that will only try to approximately solve a problem. (Yes, you could formulate the problem in terms of square of the error, then differentiate, and search for a root of the differentiated function. But that becomes a completely differeent problem. And you will still not insure a solution is found to within a tolerance. As well, you need to then eliminate any spurious solutions that may arise as maxima, or even possibly as flat spots at inflection points.)

[Mustafa]

Thank you for clarifying. Your suggestions is little complex with matlab for me. I will find a new way to achieve this.

[Torsten]

You should plot the function Ptheory1-Pexp(i,1) for a value of Pexp(i,1) for which you didn't get a solution and see if the function really becomes 0 at a certain value of the symbolic variable it depends on.

[Walter Roberson]

You could post F(106) and Vd(106) for us to test with. (Be sure to post with full significant digits)

[Mustafa]

This is F(from 100th cell to 110th cell)

1.56400000000000

1.56400000000000

1.62700000000000

1.64400000000000

1.65700000000000

1.66600000000000

1.74100000000000 (106th cell)

1.70200000000000

1.73600000000000

1.82600000000000

1.79400000000000

This is Vd (from 100th cell to 110th cell)

4303500

4342060

4446180

4419190

4484740

4527160

4546440 (106th cell)

4542580

4642850

4642850

4689120

[Walter Roberson]

Your Ptheory1 does not depend upon F or Vd. Instead, the product of F and Vd combine with other constants to become Pexp, the value that Ptheory1 must equal.

If you fplot(Ptheory1) in the area near 700, such as 650 to 750, then you will see that it peaks slightly less than your Pexp values. The peak is at approximately Te = 714.80196643093804698707623194146 with a value of approximately 1.2109926630804e-09

[Mustafa]

Thank you all for your valuable help.