I am having problem with the code written below. It is a code for finding stored strain in a polymer. Some of the errors that are coming are index exceeds matrix dimension, matrix dimensions must agree, etc. Kindly help me out.

syms epre et es(T); %stored strain

Tl = 273;
Th = 358;
Tg = 343;
Ei = 813E6; %modulus corresponding to the internal energetic deformation
N = 5.938E26; 
K = 1.38E-23; %[J/K] Boltzmann's constant
a = 2.76E-5;
n = 4;

es = zeros(1,100);

epre = .091; %pre-strain under tension

Phif = 1-1./(1+a*(Th-T).^n);

y = diff(Phif,T);

fun = @(T)1.42E-6.*T-3.16E-4;

for T = Tl:Th 
   
    Ee = 3*N*K*T;
    E = 1./((Phif./Ei)+((1-Phif)./Ee));

    et = integral(fun,Th,T); 
    
    Des = diff(es);    
    ode = Des == (epre-es-et)*y.*E./Ee;%writing stored strain as O.D.E
    cond = es(Th) == 0; %initial condition for O.D.E
    es = dsolve(ode,cond); %solving O.D.E
    
    disp(es);
    
    plot(T/343,es,'.'); hold on

end

[Geoff Hayes]

Mariam - please copy and paste the full error message to this question.

[Mariam Maroof]

Error in my coding:

Index exceeds matrix dimensions.

Error in stored_strain (line 35)

cond = es(Th) == 0; %initial condition for O.D.E

Error in run (line 96)

evalin('caller', [script ';']);

Answer 1

[Star Strider]

Try this:

syms epre et es(T); %stored strain
Tl = 273;
Th = 358;
Tg = 343;
Ei = 813E6; %modulus corresponding to the internal energetic deformation
N = 5.938E26; 
K = 1.38E-23; %[J/K] Boltzmann's constant
a = 2.76E-5;
n = 4;
epre = .091; %pre-strain under tension
    Phif = 1-1./(1+a*(Th-T).^n);
    
    y = diff(Phif,T);
    
    fun(T) = 1.42E-6.*T-3.16E-4;
   
    Ee = 3*N*K*T;
    E = 1./((Phif./Ei)+((1-Phif)./Ee));
    et = int(fun,Th,T); 
    
    Des = diff(es);    
    ode = Des == (epre-es-et)*y.*E./Ee;%writing stored strain as O.D.E
    cond = es(Th) == 0; %initial condition for O.D.E

    [VF,Sbs] = odeToVectorField(ode)
    esfcn = matlabFunction(VF)
    
    Tv = linspace(Tl, Th);
    [T,Es] = ode45(esfcn, Tv, 0);
    
    
    figure
    plot(Tv, Es)
    grid
    xlabel('T (°K)')
    ylabel('Es')

This does not give you an analytical solution (that may not actually exist), however it does give you a plot of your result.

[Mariam Maroof]

Its showing this error message:

Error using sym/subsindex (line 766)

Invalid indexing or function definition. When defining a function, ensure that the arguments are symbolic variables and the body of the function is a SYM

expression. When indexing, the input must be numeric, logical, or ':'.

Error in new (line 15)

fun(T) = 1.42E-6.*T-3.16E-4;

Error in run (line 96)

evalin('caller', [script ';']);

[Star Strider]

The code I posted ran for me without error (in R2019a).

It produced this plot:

Please check to be certain that you are running only the code I posted, and none of your previous code. (The easiest way would be to comment-out your previous code, then run only the code I posted.)

[Mariam Maroof]

The error is solved when I wrote this:

fun = @(T)1.42E-6.*T-3.16E-4;

Instead of:

fun(T) = 1.42E-6.*T-3.16E-4;

I got the plot now. Thank You.

[Walter Roberson]

You had a conflict between T being the independent symbolic variable of syms es(T), versus T being assigned a numeric value in your assignment of the output of ode45.

[Star Strider]

Mariam Maroof — As always, my pleasure. (That was the reason I suggested that you run only the code I posted.)

Walter — Thank you!