how to save only the final output of ode 45 solve
2 views (last 30 days)
Show older comments
%% 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', 1e-35, '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
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);
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:3%50%length(time_span_interval)
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),...
[0 (t_f-t_i)/2 (t_f-t_i)],initial_conditions,options);
initial_conditions=y(end,:);
end
0 Comments
Answers (1)
Sam Chak
on 23 Apr 2024
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
0 Comments
See Also
Categories
Find more on Stress and Strain 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!