I need to solve an equation involving squares and cubes.

I need to sole this equation for L, and I have a set of values for d and p.

p = -2(d/L)^3 + 3(d/L)^2

Im wondering how to put this into Matlab. Any help would be apprecitated.

[James Tursa]

Multiply both sides by L^3 and then you have a polynomial in L to work with.

Answer 1

[Dimitris Kalogiros]

It is a 3rd order equation, it has only an analytical solution.

You can use the following code:

clc; clearvars; 
 
 syms p d L
 
 p=1;
 d=3;
 
 eq= p == -2*(d/L)^3 + 3*(d/L)^2
 
 % analytical solution
 solve(eq,L)
 
 %numerical solution
 vpasolve(eq)

And in case you want to run this for a set of values , you can do the following:

 clc; clearvars;

syms p d L

% values for p, d
pValues=50:10:90;
dValues=300:50:600;

% an empty buffer for the solutions, combined with p,d
SolBuffer=cell(length(pValues)*length(dValues),3);

% solution counter
k=0;

for n=1:length(pValues)
    for m=1:length(dValues)      
        k=k+1;
        
        p=pValues(n);
        SolBuffer{k,1}=p;
        
        d=dValues(m);
        SolBuffer{k,2}=d;
        
        eq= p == -2*(d/L)^3 + 3*(d/L)^2;
        
        % analytical solution
        solve(eq, L);
        
        %numerical solution
        SolBuffer{k,3}=vpasolve(eq).';
        
    end
end


% An example: the third solution
disp('data p and d :')
disp(SolBuffer{3,1})
disp(SolBuffer{3,2})
disp('solutions:')
disp(SolBuffer{3,3})

[Dimitris Kalogiros]

Maybe, you should have to look mathworks help for more information on function root()

[Thomas McCormack]

Thank you very much, this works. Would it be easy to ask matlab to work out each combination of two sets of numbers? This is an equation for the relationship between the porosity (p), pore size (d) and the length of repeating units (L) in tissue scaffold, I am trying to work out the L values for p between 50% and 90% at 10% intervals and d values between 300 and 600 microns at 50 micron intervals in order to model them on solidworks.

I appreciate your help.

Tom

[Dimitris Kalogiros]

...I added exactly what you need... you can "accept" the answer

[Thomas McCormack]

Thats very kind of you, thanks a lot!! One final thing, is there a way to print all the answers in a matrix?

[Dimitris Kalogiros]

clc; clearvars;

syms p d L

% values for p, d
pValues=50:10:90;
dValues=300:50:600;

% an empty buffer for the solutions, combined with p,d
SolBuffer=nan(length(pValues)*length(dValues),3);

% solution counter
k=0;

for n=1:length(pValues)
    for m=1:length(dValues)      
        k=k+1;
        
        p=pValues(n);
       
        d=dValues(m);
               
        eq= p == -2*(d/L)^3 + 3*(d/L)^2;
        
        % analytical solution
        solve(eq, L);
        
        %numerical solution
        SolBuffer(k,1:3)=vpasolve(eq).';
        
    end
end


% display solutions
SolBuffer

Answer 2

[John D'Errico]

You have many values for both d and p? Do you want to solve it for all combinations of d and p?

Anyway, first, I would recognize that d and L are alwways combined in the same form, thus, we have d/L always. So first, substitute

X = d/L

Of course, once we know the value of X, we can always recover L, as:

L = d/X

Now your problem reduces to

p = -2*X^3 + 3*X^2

You could use a loop over all values of p, getting three roots for each value of p. Before we do even that however, do a plot.

ALWAYS PLOT EVERYTHING.

fun = @(X) -2*X.^3 + 3*X.^2;
fplot(fun,[-1,2])
yline(0);
yline(1);

See that I used the .^ operator there to avoid problems.

Because this is a cubic polynomial in X, you can think of your problem as having 3 real roots for X, whenever p is betweeen 0 and 1. For p larger than 1 or smaller than 0, the problem will have one real root and two complex roots.

For example,

polycoef = @(p) [-2 3 0 -p];

roots(polycoef(0))
ans =
            0
            0
          1.5
          
roots(polycoef(-.00001))
ans =
          1.5 +          0i
  -1.1111e-06 +  0.0018257i
  -1.1111e-06 -  0.0018257i

So if p is just slightly less than zero as I predicted, we see a real root, but then always two complex conjugate complex roots. The same thing happens at p==1, where 3 real roots turn into a real and two complex roots as we go above p=1.

For p between 0 and 1 however, three real roots for X.

roots(polycoef(0.5))
ans =
        1.366
          0.5
     -0.36603

Again, if you then know the value of d, we can easily recover the value of L, as

L = d./X