Equation solution is not found for all array element

2023 3 31
1 个回答
回答已采纳
更新时间：2023 4 3
8 次查看（30 天）

0 个投票

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.

请先登录，再回答此问题。

 采纳的回答

John D'Errico
John D'Errico 2023-3-31
编辑：John D'Errico 2023-3-31

0 个投票

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

10 个评论
显示 8更早的评论 隐藏 8更早的评论

Mustafa
Mustafa 2023-3-31
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
Torsten 2023-3-31
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
Mustafa 2023-3-31
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
John D'Errico 2023-3-31
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
Mustafa 2023-4-1
Thank you for clarifying. Your suggestions is little complex with matlab for me. I will find a new way to achieve this.
Torsten
Torsten 2023-4-1
编辑：Torsten 2023-4-1
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
Walter Roberson 2023-4-1
You could post F(106) and Vd(106) for us to test with. (Be sure to post with full significant digits)
Mustafa
Mustafa 2023-4-1
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
Walter Roberson 2023-4-1
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
Mustafa 2023-4-3
Thank you all for your valuable help.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Bessel functions 的更多信息

产品

版本

R2023a

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by