how do i create a for loop using symbolic matlab
3 views (last 30 days)
Show older comments
i have two symbolic variables (s ,t). Similar operations (like findinf\g the deformation gradient, euler-lagrangian strain etc) are performed using both the symbols separetly. also s varies btween 0 to t. i am creating a for loop for 's' variable where it varies between 0 to t. but i am now getting this error:
Error using :
Unable to compute number of steps from 1 to t by 1.
Error in new_14_two_syms_s_t_symbolic (line 95)
for count = 1:t
this is my code for reference:
%% SYMBOLIC
%% New method to calculate psi_r
%% scenario-1: timestep(dt) is further split into smaller time step(ds)
%% Phase 1 : Pure Deformation
%% Variable Defination
clc
clear all
close all
tic % to keep track of the time taken to run the code
%*********************** *********************
%% Constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta = 80+273; % absloute temperature in Kelvin
c_small = 1; %homogenous oxygen penetration
R = 8.314; % universal gas constant in J/mol
%***********************************************
nu_r = 8.745e4; % in /sec
E_r = 9.5852e4; % in J/mol
cons_r = c_small * nu_r * exp(-(E_r / (R * theta)));
% %*******************************************************
%% SYMBOLIC VARIABLES DEFINATION
syms t s
% %*******************************************************
%% Variable Defination
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % for an aged specimen
lam_0 = 1; %initial value of stretch , i.e stretch at t=0
dlam= 0.25*10^(-3); %dlam is rate of stretch %lam_t=dlam*t+lam_0;
% s=0;% value of "smaller timestep". this value gets updated at the beginning of each loop (counter variable is 'i_s').
%% Evalaution of quantities BEGINS HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lam_t = dlam * t + lam_0; %lam_t=at+b , at t=0, lam_t=1 (i.e strain=0%)
strain_t=(lam_t-1)*100; %strain percentage
% q_r(i) = 1 - exp(-cons_r * t);
%% Tensoral quantities
% F-->Deformation gradient ,lam is short for lamda or the stretch
% F = zeros(3, 3);
for ii = 1:3
for jj = 1:3
if ii == jj
if ii == 2
F(ii, jj) = (dlam)*t + lam_0;
elseif ii == 1 || ii == 3
F(ii, jj) = 1/((dlam*t + lam_0)^0.5);
end
end
end
end
%********************************************
% Fdot-->is the differentiation of Deformation gradient ,lam is short for lamda or the stretch
% Fdot = zeros(3, 3);
for m = 1:3
for n = 1:3
if m == n
if m == 2
Fdot(m, n) = dlam;
elseif m == 1 || m == 3
Fdot(m, n) = -(0.5*dlam)/((dlam*t + lam_0)^1.5);
end
end
end
end
%********************************************
% C -->Cauchy-Green deformation tensor
C = transpose(F)*F;
%********************************************
% Cdot -->is the differentiation of Cauchy-Green deformation tensor
Cdot = transpose(Fdot)*F + Fdot*transpose(F);
%********************************************
% E -->Green-Lagrangian Strain
E = 0.5*(C - eye(3));
%********************************************
J = 1; % Isochoric
mu0 = 1.844e6;
K = mu0*1e3; % Bulk Modulus
%********************************************
C_inv = inv(C);
Ic = trace(C);
I = eye(3);
mu_r = 567.3e6;
% mu_r = 5.67e6;
%********************************************
%********************************************
for count = 1:t
% parameters for inner loop
%********************************************
dlam_s=dlam;% value by which stretch (dlam_st) increases for every iteration of the loop
lam_s = dlam_s * s + lam_0;
s_num=10;%s_num is the number of "smaller timesteps"
ds=0.1;%ds is size of each "smaller timestep"
term_A2= zeros(3,3);
term_A= zeros(3,3);
% term_TOTAL = zeros(1,time_steps); %COMMENTED only in this symbolic version not in the numerical one.
%*******************************************
% F_s = zeros(3, 3);
for ii_s = 1:3
for jj_s = 1:3
if ii_s == jj_s
if ii_s == 2
F_s(ii_s, jj_s) = (dlam_s)*s + lam_0;% F is same in both the outer and inner loops . RIGHT ?? CONFIRM THIS !!!
elseif ii_s == 1 || ii_s == 3
F_s(ii_s, jj_s) = 1/((dlam_s*s + lam_0)^0.5);
end
end
end
end
%********************************************
% F_sdot-->is the differentiation of Deformation gradient ,lam is short for lamda or the stretch
% F_sdot = zeros(3, 3);
for m_s = 1:3
for n_s = 1:3
if m_s == n_s
if m_s == 2
F_sdot(m_s, n_s) = dlam_s;
elseif m_s == 1 || m_s == 3
F_sdot(m_s, n_s) = -(0.5*dlam_s)/((dlam_s*s + lam_0)^1.5);% Fdot is same in both the outer and inner loops . RIGHT ?? CONFIRM THIS !!!
end
end
end
end
%********************************************
% C -->Cauchy-Green deformation tensor
C_s = transpose(F_s)*F_s;
%********************************************
% Cdot -->is the differentiation of Cauchy-Green deformation tensor
C_sdot = transpose(F_sdot)*F_s + F_sdot*transpose(F_s);
%********************************************
% E -->Green-Lagrangian Strain
E_s = 0.5*(C_s - eye(3));
%********************************************
J = 1; % Isochoric
mu0 = 1.844e6;
K = mu0*1e3; % Bulk Modulus
%********************************************
C_s_inv = inv(C_s);
Ic_s = trace(C_s);
I = eye(3);
mu_r = 567.3e6;
% mu_r = 5.67e6;
%D1_tp = tensorprod(1/3 * C_inv , C_inv); %D1 obtained using tensorprod function
%D1 = zeros(3, 3, 3, 3); % should i initialize D1 as an identity tensor of 3x3x3x3 first?
for i1 = 1:3
for j1 = 1:3
for k1 = 1:3
for l1 = 1:3
D1(i1,j1,k1,l1) = (1/3) * C_inv(i1,j1) * C_inv(k1,l1);
end
end
end
end
%D2_tp = tensorprod(C_inv , C_inv); %D2 obtained using tensorprod function
%D2 = zeros(3, 3, 3, 3);
for i2 = 1:3
for j2 = 1:3
for k2 = 1:3
for l2 = 1:3
D2(i2,j2,k2,l2) = C_inv(i2,k2) * C_inv(j2,l2); %Notice the indices (i,k) and (j,l). The transpose is done between 2nd index of first C_inv and 1st index of second C_inv (T23 as said in paper).This straight away ensures the special transposition in lines below
end
end
end
end
%Special transposition done below------DOUBT
% a=D2_tp(2,3,:,:);
% b=D2_tp(3,2,:,:);
% D2_tp(3,2,:,:)=a;
% D2_tp(2,3,:,:)=b;
%D3_tp = tensorprod(-1*C_inv , I);%D3 obtained using tensorprod function
%D3 = zeros(3, 3, 3, 3);
for i3 = 1:3
for j3 = 1:3
for k3 = 1:3
for l3 = 1:3
D3(i3,j3,k3,l3) = (-1) * C_inv(i3,j3) * I(k3,l3);
end
end
end
end
%D4_tp = tensorprod(-I , C_inv); %D4 obtained using tensorprod function
%D4 = zeros(3, 3, 3, 3);
for i4 = 1:3
for j4 = 1:3
for k4 = 1:3
for l4 = 1:3
D4(i4,j4,k4,l4) = (-1)* I(i4,j4) * C_inv(k4,l4);
end
end
end
end
cons = 1/6 * mu_r * J^(-2/3);
dd_sed = cons * (Ic * (D1 + D2) + D3 + D4);
D = 4 * dd_sed;
q_r_s = 1 - exp(-cons_r * s);
term_x=D.*q_r_s;%term_x is an intermediate term.
%term_A1=transpose(term_x);
% Perform the transpose(t) using nested loops
for i_tran = 1:3
for j_tran = 1:3
for k_tran = 1:3
for l_tran = 1:3
term_A1(i_tran, j_tran, l_tran, k_tran) = term_x(i_tran, j_tran, k_tran, l_tran);% only the last two indices, seem to transpose. PLZ CONFIRM THIS !!!
end
end
end
end
term_A2_interim=E-E_s; % should i take E(i), i.e at a particular value of outer loop iteration (i) or simply E ??
term_A2 = term_A2 + term_A2_interim;
% term_A = zeros(3,3);
% Compute the tensor contraction between term_A1 and term_A2 to obtain term_A
%iA,jA,kA,lA are loop counter variables involved in this 'Contraction'
%process.
for iA = 1:3
for jA = 1:3
for kA = 1:3
for lA = 1:3
term_A_interim(iA,jA) = term_A1(iA,jA,kA,lA) * term_A2(kA,lA);%the last two indices of A1 and A2 are contracted with each other
% IMP---this expression is different from the corresponding numerical matlab version file
end
end
end
end
term_A=term_A+term_A_interim;
term_B = term_A2;
% term_TOTAL(i_s) = 0;% at the end, term_TOTAL will be a scalar , so should i make a blank 3*3 matrix for it?
% Compute the tensor contraction between term_A and term_B to obtain term_TOTAL
%iTOT,jTOT,kTOT,lTOT are loop counter variables involved in this 'Contraction'
%process.
% for iTOT = 1:3
% for jTOT = 1:3
% term_TOTAL(i_s) = (term_A(iTOT,jTOT) * term_B(iTOT,jTOT));%the last two indices of D and Cdot are contracted with each other
% % at the end, term_TOTAL will be a scalar , so should i simply write term_TOTAL instead of term_TOTAL(iTOT,jTOT)
% % IMP---this expression is different from the corresponding numerical matlab version file
% end
% end
term_TOTAL_interim=(term_A.*term_B);
term_TOTAL=sum(term_TOTAL_interim(:));
end
% %% INNER LOOP ENDS HERE
A=term_TOTAL;
A_after_int=int(A,s);
A_after_limits=subs(A_after_int,s,t)-subs(A_after_int,s,0);
%% Evaluation of quantities ENDS HERE
%% Plots
lw=8; % linewidth
ms=25; % marker size
dp=10^4; % marker every 'dp' Data Point . it is calculated as (0.1*time_steps) , such that we can get exactly 10 data points
fs= 50; % font size
fsl=30; % font size in Legend
%S all components along X-----------------
figure()
xlim([0 10]);
%% Save *.mat file
% %********************************************
save("Case_Fig_1_Y_10_7_Sec_80degC_phase_1_SYMBOLIC.mat")
% save("Case_Fig_1_Y_10_7_Sec_80degC.csv")
elapsedTime = toc;% to keep track of the time elapsed to run the code
% Display the elapsed time
disp(['Elapsed time: ' num2str(elapsedTime) ' seconds']);
3 Comments
Steven Lord
on 17 Jan 2024
Lets say t varies from 0 to 1000 seconds, then at say t=700th sec, 's' varies from 0 to 700. similarly at say t= 800th seconds varies from 0 to 800 and so on.
So don't use the symbolic variable t in your for loop. Use a numeric value as the limit of your loop, possibly with subs to substitute the numeric value for t in your expression.
Answers (0)
See Also
Categories
Find more on Calculus 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!