I created a dummy model to test the for-loop code.
Though you didn't annotate the code properly, I believe that time span should be
[t_i, t_f]
instead of
[0 (t_f-t_i)/2 (t_f-t_i)]
%% RSFODE-CRACK_III QUASI-STATIC APPROXIMATION
% Parameter values as in Rubin & Ampuero (2005)
tic
V_init=10e-9; %in meters/s
Tau_dot_inf=10e-2; % Pa/s remote stressing rate
sigma=100; %Mpa normal stress
mu_dash=11.56; %Gpa rigidity modulus.
options = odeset('RelTol', 3e-14, 'AbsTol', 1e-6, 'MaxStep', 1e5);
Gmod = 3e4; % in MPa
v_shear = 3e9; % in microns/sec
eta= 0.5*Gmod/v_shear; % Radiation damping term
%RATE AND STATE PARAMETERS
D_c=1; %in mm
b=0.01; % rate state parameter that characterize the magnitude of evolution effect.
a_b=0.5; % a/b ratio
a=a_b*b; % Rate state parameter that characterize the magnitude of direct effect.
%FAULT LENGTH AND DISCRETIZATION ELEMENT
% Nx=2^12; % No. of discretization element
Nx=3;
spacing=0.05; % spacing between the points
Lx=Nx*spacing; % Length of the fault in x-direction.
% Creating a random distribution of v values
V_init_dist=10e-12 + (10e-11-10e-12).*rand(Nx,1);
%kernel
Kx = 2*pi/Nx * [0:Nx/2-1, 0, -Nx/2+1:-1];
%parameters that need to be passed on to the derivate script
par = [a,b,D_c,sigma,Tau_dot_inf,mu_dash];
law = 'A';
% create time intervals for the iteration
% time_span_interval=linspace(0,1e9,2e7);
time_span_interval=linspace(0, 60, 61);
t_initial=0; % Initial time
% Create a distribution of theta
theta_init=(D_c./V_init_dist);
% PREALLOCATING a single vector for initial conditions
initial_conditions = zeros(2*Nx, 1);
initial_conditions(1:Nx,1) = theta_init;
initial_conditions(Nx+1:2*Nx,1) = V_init_dist;
% SOLVING FOR V AND THETA
for i = 1:length(time_span_interval)-1
t_i = time_span_interval(i);
t_f = time_span_interval(i+1);
[t,y] = ode45(@(t,y) rsfode_crack_III(t, y, Nx, par, law), [t_i t_f], initial_conditions, options);
initial_conditions=y(end,:);
plot(t, y), hold on
end
% [t,y] = ode45(@(t,y) rsfode_crack_III(t,y,Nx,par,law), time_span_interval, initial_conditions, options);
% plot(t, y), hold on
grid on
%% Dummy ODE model
function dy = rsfode_crack_III(t,y,Nx,par,law)
A = [zeros(2*Nx-1, 1), eye(2*Nx-1);
- par];
B = [zeros(2*Nx-1, 1); 1];
K = lqr(A, B, eye(2*Nx), 1);
u = - K*y;
dy = A*y + B*u;
end
