Hello everyone,
I am trying to solve an ordinary differential equation that has a function in it that only activates when T=0.05 (essentially a dirac impulse function). I need o figure out a way to ensure that even with an adaptive time step this function will trigger. Below I've attached my ode function and my script to run the file. Any help would be greatly appreciated.
(Please note that I do not use the actual dirac function in the ode because it causes values to go to NaN because of the linear algebra involved. My ode function:
function [Xddot]=mdof(T,X)%System Parameters
m=30; %kg
M=200; %kga=3; %m
b=1; %mk1=400*1000; %N/m
k2=k1;k3=100*1000; %N/mk4=k3;I=200; %kg*m^2
C1=3200; %Ns/m
C2=C1;C3=800; %Ns/mC4=C3;M2=[(M*b^2+I)/(a+b)^2 (M*a*b-I)/(a+b)^2 0 0;... (M*a*b-I)/(a+b)^2 (M*a^2+I)/(a+b)^2 0 0;... 0 0 m 0;... 0 0 0 m];C=[ C1 0 -C1 0 ;... 0 C2 0 -C2;... -C1 0 C1+C3 0 ;... 0 -C2 0 C2+C4];K=[ k1 0 -k1 0 ;... 0 k2 0 -k2;... -k1 0 (k1+k3) 0 ;... 0 -k2 0 (k2+k4)]; Xdot=[X(5);X(6);X(7);X(8)]; XX=[X(1);X(2);X(3);X(4)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F=[0;0;(T==0)*1e100;((T-0.05)==0)*1e100] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T Z=(F-C*Xdot-K*XX); Xddot2=M2\Z; Xddot=[Xdot;Xddot2]; end
Here the part set off is where I am having the problem. The adaptive time step jumps over T=0.05 and thus never triggers the input. If i try to use a error tolerance though, such as
(abs(T-0.05)<=1e-2)
I either end up with the same problem of the function never triggering or a problem of the function triggering numerous times because the time step slows way down to handle the high input.
Here's the call function for if you want to run it.
clc, clear, close all%%Numerical Method
IC=zeros(1,8);t=0:0.0001:5;options=odeset('RelTol',1e-6,'AbsTol',1e-6);[T,X]=ode45('mdof',t,IC,options);figureplot(T,X(:,1),'r-')hold onplot(T,X(:,2),'k-')plot(T,X(:,3),'g-')plot(T,X(:,4),'-')title('\textbf{Numeric Analysis Results}','interpreter','LaTeX','FontSize',16)xlabel('Time, s','interpreter','LaTeX','FontSize',14)ylabel('Displacement, m','interpreter','LaTeX','FontSize',14)legend('x_1','x_2','x_3','x_4')hold off
Thanks so much for your help!!
Best Answer