Trial and error problem

Dear Matlab Community Members

My objective is to find the values a, b, and c. of the following equation:

log(m)=log(a)+(b/(T-c))

I have three values of m & T

m=0.0701 & T=293.65

m=0.0262 & T=313.15

m=0.00433 & T=373.15

I want to know what is the best code or technique to find these values.

Thank you very much

[Torsten]

So you have one data point m=0.0262, T=313.15 and you want to find parameters a, b and c such that log(m)=log(a)+(b/(T-c)) ?

[Ahmad Al-Issa]

yes, you are right.

[Torsten]

Why not a = m, b = 0, c = arbitrary ?

But seriously: If you have many data points for m and T, you can use a fitting tool to fit the constants.

But adjusting three parameters with only one data point doesn't make sense since the (almost arbitrary) parameters for this single data point will most probably not work for a second one.

[Ahmad Al-Issa]

actually, I use the value of m and T only for one time because it is measured experimentally.

after finding a, b, and c, I will use the equation to find m at any T.

could you please give me the most efficient code to execute the trial-and-error method? this is the only way to solve the problem, I guess.

[Torsten]

actually, I use the value of m and T only for one time because it is measured experimentally.

after finding a, b, and c, I will use the equation to find m at any T.

As I tried to explain in my response, exactly this doesn't make sense since you only have one measurement point.

[Ahmad Al-Issa]

yes you are right.

I change my post. there are another values for m & T. you can see the post again.

that is mean now we have 3 equation and 3 unknown variables a, b, and c.

so what is the best code to find them?

thank you in advanced

[the cyclist]

How did you choose that function for the fit? It seems like a poor one, because of the singularity that will happen when T==c.

c = 280;
T = 260 : 380;
f = @(t) 1./(t -c);

plot(T,f(T))

Even if your values of T are strictly greater than c, (the "critical" temperature), this behavior will make it difficult for any algorithm to fit it.

One can try (e.g. with fitnlm), but MATLAB spits out warnings about overparameterization, and the fit is not good.

[Ahmad Al-Issa]

@the cyclist

thank you for your answer, but this is not right.

I found the soution by solve a system of equations.

m1=0.0701;
T1=293.65;
m2=0.0262;  
T2=313.15;
m3=0.00433;  
T3=373.15;

syms a b c
eqns = [log(a)+(b/(T1-c))== log(m1), log(a)+(b/(T2-c))== log(m2), log(a)+(b/(T3-c))== log(m3)];
S = (vpasolve(eqns,[a b c]));

Answer 1

[Torsten]

The results are not convincing, but that's the way to go with three data points.

fun1 = @(x)[0.0701-x(1)*exp(x(2)/(293.65-x(3)));0.0262-x(1)*exp(x(2)/(313.15-x(3)));0.00433-x(1)*exp(x(2)/(373.15-x(3)))];
fun2 = @(x)[log(0.0701)-(log(x(1))+x(2)/(293.65-x(3)));log(0.0262)-(log(x(1))+x(2)/(313.15-x(3)));log(0.00433)-(log(x(1))+x(2)/(373.15-x(3)))];
x0 = [2 1 450];
x = fsolve(fun1,x0)

No solution found.

fsolve stopped because the problem appears regular as measured by the gradient,
but the vector of function values is not near zero as measured by the
value of the function tolerance.

x = 1×3

0.0345 8.2224 633.4040

fun1(x)

ans = 3×1

0.0364 -0.0075 -0.0291

x = fsolve(fun2,x0)

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+02.

x = 1×3

0.0217 53.9522 974.0568

fun2(x)

ans = 3×1

1.2519 0.2701 -1.5220

[the cyclist]

"Not convincing", indeed. :-)

% The data
m = [0.0701;
     0.0262;
     0.0043];

T = [293.65;
     313.15;
     373.15];

% Solve for the coefficients at the data points
x0 = [2 1 450];
fun2 = @(x)[log(m(1))-(log(x(1))+x(2)/(T(1)-x(3)));
            log(m(2))-(log(x(1))+x(2)/(T(2)-x(3)));
            log(m(3))-(log(x(1))+x(2)/(T(3)-x(3)))];
x = fsolve(fun2,x0) % these are the coefficents

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+02.

x = 1×3

0.0217 54.0870 974.0394

% The resulting fitted function to get log(m) from T
fun_fsolve = @(t) log(x(1))+x(2)./(t-x(3));

% Calculate the predicted values in the region of the data
T0 = 290:390;
logm_fsolve = fun_fsolve(T0);

% Plot the data and the fit
figure
hold on
plot(T,log(m),'*')
plot(T0,logm_fsolve,'r')
legend('data','fsolve')

[Ahmad Al-Issa]

@Torsten

thank you for your answer.

I found the solution by solving a system of equations.

m1=0.0701;
T1=293.65;
m2=0.0262;  
T2=313.15;
m3=0.00433;  
T3=373.15;
sym abc
eqns = [log(a)+(b/(T1-c))== log(m1), log(a)+(b/(T2-c)))== log(m2), log(a)+ (b/(T3-c))==log(m3)];
S = (vpasolve(eqns,[abc]));

[Torsten]

Seems to work - I'm surprised about the symbolic solver:

m1=0.0701;
T1=293.65;
m2=0.0262;  
T2=313.15;
m3=0.00433;  
T3=373.15;
syms a b c
eqns = [log(a)+b/(T1-c)== log(m1), log(a)+b/(T2-c)== log(m2), log(a)+ b/(T3-c)==log(m3)];
S = vpasolve(eqns,[a b c])

S = struct with fields:
    a: 0.00007308300028591661380690628944468
    b: 800.19428469747439653736436104991
    c: 177.10693349035606680061157547077

double(subs([log(a)+b/(T1-c)- log(m1), log(a)+b/(T2-c)- log(m2), log(a)+ b/(T3-c)-log(m3)],[a b c],[S.a,S.b,S.c]))

ans = 1×3

1.0e-38 * 0 -0.1469 -0.1469

[Ahmad Al-Issa]

@Torsten

@the cyclist

thanks for your try.

The right answer, as it is in the refrences, is:

a=0.000073

b=800

c=177

because a, b, and c are constants and they need 3 equations to solve.

[Torsten]

@AHMAD AL-ISSA

Isn't the form of the equation usually

m = 10^(A-B/(C+T))

where T is in degreeC ?

Thus

log10(m) =  A - B/(T+C)

where log10 is not the natural, but the decadic logarithm and T is temperature in degreeC ?

m1=0.0701;
T1=293.65-273.15;
m2=0.0262;  
T2=313.15-273.15;
m3=0.00433;  
T3=373.15-273.15;
syms a b c
eqns = [a-b/(T1+c)== log10(m1), a-b/(T2+c)== log10(m2), a-b/(T3+c)==log10(m3)];
S = vpasolve(eqns,[a b c])

S = struct with fields:
    a: -4.1361836318056620904194888300452
    b: -347.51996229463221164872269115385
    c: 96.04306650964382392391269916326

double(subs([a-b/(T1+c)- log10(m1), a-b/(T2+c)- log10(m2), a-b/(T3+c)-log10(m3)],[a,b,c],[S.a,S.b,S.c]))

ans = 1×3

1.0e-38 * 0.0735 0.0735 0.1469