In this system, there are several problems, including:

1- The indices for the two systems cannot be equal.

2- The global function is only valid for one value of (MD) (the last value), while I need all values from the first system.

68 views (last 30 days)

Show older comments

mohammed
on 20 May 2024 at 10:02

function [sol1] = two_system

clc

close all

clear all

global z; %%problems

%%1- The indices for the two systems cannot be equal.

%%2- The global function is only valid for one value of (MD) ...

%%(the last value), while I need all values from the first system.

z= 0.05;

tau = 0.1;

tspan = linspace(0, 5, 300);

opt = odeset('RelTol',1e-3);

sol1 = dde23(@system1,tau ,[0],tspan, opt);

[x1, y1] = ode23(@system2, tspan, [0], opt);

figure(1)

plot(sol1.x, sol1.y(1,:),'--r','LineWidth',1.5)

hold on

plot(x1, y1(:,1),'b','LineWidth',1.5)

L2=length(x1)

L1=length(sol1.x)

end

function dydt = system1(t,y,D)

global z MD;

M1 = y(1);

MD = D(1);

s1 = MD %for checking

F_M1= 2;

theta1=0*pi/180;

A1 = 0.01;

%%%%%%% equation_system1 %%%%%%%

dMdt=(A1*(2*pi*(F_M1))*cos(2*pi*(F_M1)*(t)+theta1));

dydt = [dMdt];

end

function dydt = system2(t,y)

global z MD;

s2=MD %% for checking

M2 = y(1);

F_M2= 2;

theta2=0*pi/180;

A2 = 0.01;

%%%%%%% equation_system2 %%%%%%%

dMdt=(A2*(2*pi*(F_M2))*cos(2*pi*(F_M2)*(t)+theta2))+MD*z;

dydt = [dMdt];

end

Torsten
on 20 May 2024 at 20:02

Edited: Torsten
on 20 May 2024 at 20:49

I can't follow your description. Maybe I understand better if you write down the equations for the two systems you try to solve in a mathematical way.

Which value for MD from the solution of the first system do you want to insert for the second system in the line

dMdt=(A2*(2*pi*(F_M2))*cos(2*pi*(F_M2)*(t)+theta2))+MD*z;

?

If it is the solution of the first equation, evaluated at t-tau, I don't see the problem solving both equations simultaneously. If you don't use this term in the definition of first equation, the output of the first system won't change with tau, either.

@Sam Chak did exactly this in the code he suggested.

Sam Chak
on 20 May 2024 at 16:44

Hi @mohammed

Based on the information you have provided, it appears the two systems could potentially be combined and then solved using the dde23 solver. However, I am uncertain about the specific approach you would like to use to calculate the correlation between and .

%% call dde23 solver

tau = 0.1;

tspan = [0 1];

sol = dde23(@ddefun, tau, @history, tspan);

t = sol.x;

M = sol.y;

%% Check the lengths of solutions M1 and M2

size(M)

%% Plot results

subplot(211)

plot(t, M(1,:), '-o', 'linewidth', 1.5, 'color', '#F09494'), grid on

xlabel('Time t'), ylabel('Amplitude'), title('Solution M1')

subplot(212)

plot(t, M(2,:), '-o', 'linewidth', 1.5, 'color', '#AAE2A8'), grid on

xlabel('Time t'), ylabel('Amplitude'), title('Solution M2')

%% Two-system DDE

function dMdt = ddefun(t, M, Z)

M1lag = Z(:,1); % delayed state M1

F_M1 = 2;

theta1 = 0*pi/180;

A1 = 0.01;

F_M2 = 2;

theta2 = 0*pi/180;

A2 = 0.01;

z = 0.05;

dMdt = [A1*2*pi*F_M1*cos(2*pi*F_M1*t + theta1);

A2*2*pi*F_M2*cos(2*pi*F_M2*t + theta2) + M1lag(1)*z];

end

%% historical data

function s = history(t)

s = [0;

0];

end

Sam Chak
on 21 May 2024 at 12:02

Hi @mohammed

Both @Torsten and I have been trying to understand the problem you are attempting to solve. However, it appears the selected value for the factor z is too small, causing the product term 'M1lag(1)*z' to become insignificant in your toy problem.

Would it be possible to instead run a simulation on the chaotic van der Pol oscillator, a double system? I can provide an initial setup with a constant delay of 20 seconds and a factor of 1. You can then adjust the parameters as needed and explain the specific observations you would like to make.

Exploring the chaotic van der Pol oscillator system may provide more meaningful insights than the current approach.

%% call dde23 solver

tau = 20; % time delay

tspan = [0 40];

sol = dde23(@ddefun, tau, @history, tspan);

t = sol.x;

x = sol.y;

%% Plot results

tL = tiledlayout(2, 2, 'TileSpacing', 'Compact');

nexttile

plot(t, x(1,:), 'linewidth', 1.5, 'color', '#115AD0'), grid on, ylim([-3 3])

xlabel t, ylabel Amplitude, title('Response, x1(t)')

nexttile

plot(t, x(3,:), 'linewidth', 1.5, 'color', '#EC3114'), grid on, ylim([-3 3])

xlabel t, ylabel Amplitude, title('Response, x3(t)')

nexttile

plot(x(1,:), x(2,:), 'linewidth', 1.5, 'color', '#F1BC5F'), grid on, axis([-3 3 -5 5])

xlabel x_{1}, ylabel x_{2}

nexttile

plot(x(3,:), x(4,:), 'linewidth', 1.5, 'color', '#3A6A1E'), grid on, axis([-3 3 -5 5])

xlabel x_{3}, ylabel x_{4}

%% Double van der Pol oscillator

function dxdt = ddefun(t, x, Z)

x1Lag = Z(:,1); % delayed state x1

factor = 1; % similar to z in your System 2

% ODEs

dxdt = [x(2);

- x(1) + (1 - x(1)^2)*x(2);

x(4);

- x(3) + (1 - x(3)^2)*x(4) + x1Lag(1)*factor];

end

%% historical data

function s = history(t)

s = [1; % initial value

0;

1;

0];

end

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

Start Hunting!