how do i create a for loop using symbolic matlab
3 次查看(过去 30 天)
显示 更早的评论
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 个评论
Steven Lord
2024-1-17
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.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!