find unknow value to plot(best fit)

1 view (last 30 days)
Carc
Carc on 18 Jan 2023
Commented: Walter Roberson on 19 Jan 2023
I would like to find unknown values(D a b) below the equation.
Could you let me know what errors my code have. Finally, I am going to find the best fit curve. Thank you.
data = readtable('OO_band');
x = data.energy; %% known from data
r = data.distance; %% known from data
syms D a b; %% unknown
potential = x == D*((1-Exp(-a*(r-b)))^2 - 1);
sol = solve(potential, [D a b]);

Answers (1)

Walter Roberson
Walter Roberson on 19 Jan 2023
Your x and r are vectors, so your potential is going to create a vector of equations. You are then trying to solve() that vector of equations for values of D, a, b that exactly satisfy all of the equations simultaneously . That is unlikely to work unless you happen to have exactly 3 rows in your file.
For best fit you are not looking for coefficients that solve() all of the equations exactly, you are looking for coefficients that give the least overall error, for some meaning of error . For example you might look for smallest sum of squares of error.
rng(123456)
N = 4;
x = rand(N,1) * 10;
r = rand(N,1) * 10;
syms D a b; %% unknown
model = D*((1-exp(-a*(r-b))).^2 - 1)
model = 
actual = x
actual = 4×1
1.2697 9.6672 2.6048 8.9724
residue = expand(sum((model - actual).^2))
residue = 
syms AB
residueAB = subs(residue, a, AB/b)
residueAB = 
partial_D = solve(diff(residueAB, D),D)
partial_D = 
residue2 = subs(residueAB, D, partial_D)
residue2 = 
symvar(residue2)
ans = 
db = diff(residue2,b)
db = 
partial_b = solve(db,b)
Warning: Unable to find explicit solution. For options, see help.
partial_b = Empty sym: 0-by-1
You can see from this that pure symbolic solutions are not going to work, so ...
partial = vpasolve(db)
partial = struct with fields:
AB: -249.80509071384312325902061445948 b: 10.187373506758686423500301869879
subs(partial_D, partial)
ans = 
  1 Comment
Walter Roberson
Walter Roberson on 19 Jan 2023
After which you would need to convert AB and B back to A
(The equations had exp(constant*a*b) in them a lot; I was hoping that by combining them that it might make it easier to solve for the variable, but it did not, so probably you could work more directly on the original equation after substituting in the partial solution for D.)

Sign in to comment.

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!