If I have a sigmoid equation depend on time and I want to input to Runge Kutta Orde 4Th equation. So, the graph show the sigmoid from Runge_kutta . How to call the function?

1 view (last 30 days)
This is my sigmoid equation depend on time ( j is looping for time)
%input for time
t(1)=0;
dt=0.1; %time interval
t=0:dt:100; %time span
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
for j = 1:length(t)-1
T(j+1)=T(j)+dt;
M1(j+1)= M1(j)+1./(1+exp(-T(j))); %this is sigmoid equation
end
figure
plot (T,M1,'r','Linewidth',3)
xlabel('time')
ylabel('M1')
That equation I want to input into Runge kutte orde 4 and show the sigmoid graph from RK 4
Thank you

Accepted Answer

Sam Chak
Sam Chak on 29 May 2023
Edited: Sam Chak on 29 May 2023
I'm unsure if I understand the issue correctly. Do you expect a plot like this?
tspan = linspace(0, 100, 10001);
M0 = 0; % initial value
[t, M] = ode45(@odefcn, tspan, M0);
% Solution Plot
plot(t, M, 'linewidth', 1.5),
grid on, xlabel('t'), ylabel('M(t)')
% System
function Mdot = odefcn(t, M)
Mdot = 1/(1 + exp(- t));
end
  2 Comments
Sam Chak
Sam Chak on 29 May 2023
The integration method in your code looks like the forward Euler method and it produce similar result like the RK4.
t(1) = 0;
dt = 1e-2; % time interval
t = 0:dt:100; % time span
T = zeros(length(t) + 1, 1); % empty array for t
M1 = zeros(length(t) + 1, 1); % empty array for M1
for j = 1:length(t)-1
T(j+1) = T(j) + dt;
M1(j+1) = M1(j) + dt./(1 + exp(- T(j)));
end
figure
plot(t, M1(1:end-1), '--'), grid on
xlabel('time')
ylabel('M1')

Sign in to comment.

More Answers (1)

Sam Chak
Sam Chak on 30 May 2023
The RK4 algorithm by Cleve Moler, is available in File Exchange.
t = linspace(0, 100, 1001);
M0 = 0; % initial value
Mout = ode4(@odefcn, 0, 0.1, 100, M0);
% Solution Plot
plot(t, Mout, 'linewidth', 1.5),
grid on, xlabel('t'), ylabel('M(t)')
% System
function Mdot = odefcn(t, M)
Mdot = 1/(1 + exp(- t));
end
% Classical Runge-Kutta 4th-order ODE solver
function yout = ode4(F, t0, h, tfinal, y0)
y = y0;
yout = y;
for t = t0 : h : tfinal-h
s1 = F(t,y);
s2 = F(t+h/2, y+h*s1/2);
s3 = F(t+h/2, y+h*s2/2);
s4 = F(t+h, y+h*s3);
y = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
yout = [yout; y];
end
end

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!