how do i create a for loop using symbolic matlab
    5 visualizzazioni (ultimi 30 giorni)
  
       Mostra commenti meno recenti
    
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 Commenti
  Steven Lord
    
      
 il 17 Gen 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.
Risposte (0)
Vedere anche
Categorie
				Scopri di più su Calculus in Help Center e File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!


