I am trying to solve this coupled ODE from the past few days using ode45 but there's an error at ode45 function which I am unable to understand .

Question

Yogesh 2024-1-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2065876-i-am-trying-to-solve-this-coupled-ode-from-the-past-few-days-using-ode45-but-there-s-an-error-at-ode

评论： Yogesh 2024-1-3

采纳的回答： Sam Chak

在 MATLAB Online 中打开

close all
clear all
clc
alpha=3.55*10^-3;
sigma=3.7292*10^11;
eta=rand;
TB=3.59*10^8;
delta_t=50*10^-12;
n=1.45;
c=3*10^8;
Q=141.3227*10^-4;
syms t p(t) AL(t) AS(t) f Y
Eq1=diff(AL)==1i*sigma*p(t)*AS(t);
Eq2=diff(AS)==1i*sigma*conj(p(t))*AL(t);
Eq3=(alpha/TB)*diff(p,2)+(alpha-1i)*diff(p)-((1i*TB)/2)*p(t)==eta*AL(t)*conj(AS(t))-1i*f;
[VF,subs]=odeToVectorField(Eq1,Eq2,Eq3);
disp(VF)
% disp(subs)
ftotal=matlabFunction(VF);
real_part=randn();
imag_part=randn();
complex=exp(-(real_part+1i*imag_part))^2;
f=sqrt((n*Q)/(delta_t)^2*c)*complex;
r=sqrt((n*Q)/(c*TB))*complex;
ic=[0 0 r];
tspan=[1 2];
[t, sol] = ode45(@(t,Y) ftotal(t,Y), tspan, ic);
plot(t, sol(:, 1), 'LineWidth', 2, 'DisplayName', 'Eq1');
hold on;
plot(t, sol(:, 2), 'LineWidth', 2, 'DisplayName', 'Eq2');
plot(t, sol(:, 3), 'LineWidth', 2, 'DisplayName', 'Eq3');
xlabel('Time');
ylabel('Solution');
title('Solution of the Coupled ODE System');
grid on;

I have attached the question and the code which I am using it to solve the equations but there is an error at ode45 which I cant understand.

2 个评论
显示无隐藏无

Maven Euler 2024-1-3

could you please paste your code in the code snippet? It is cumbersome to do it from the png.

Sulaymon Eshkabilov 2024-1-3

Please post your whole code. The community can simulate your code and correct it instead of re-typing the whole code.

It looks like you have got a heavy part done.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2024-1-3

在 MATLAB Online 中打开

Hi @Yogesh

The first issue is that the free parameter "f" must be defined before its inclusion in the ODEs. All constants and parameters should be moved before the ODE section. The second issue is that since there are 4 states, 4 initial values are required, but only 3 initial values are provided. Theoretically, the code should run after fixing these two main issues. However, the ODE solver fails to integrate from the beginning. This is likely caused by the complex-valued ODEs. Please check the equations.

%% Constant

alpha=3.55*10^-3;

sigma=3.7292*10^11;

eta=rand;

TB=3.59*10^8;

delta_t=50*10^-12;

n=1.45;

c=3*10^8;

Q=141.3227*10^-4;

syms t p(t) AL(t) AS(t) f Y

% --- move to here ---

%% Parameters

real_part=randn();

imag_part=randn();

complex=exp(-(real_part+1i*imag_part))^2;

f=sqrt((n*Q)/(delta_t)^2*c)*complex; % <-- this parameter must be defined before ODE

r=sqrt((n*Q)/(c*TB))*complex;

% --- move to here ---

%% ODE section

Eq1=diff(AL)==1i*sigma*p(t)*AS(t);

Eq2=diff(AS)==1i*sigma*conj(p(t))*AL(t);

Eq3=(alpha/TB)*diff(p,2)+(alpha-1i)*diff(p)-((1i*TB)/2)*p(t)==eta*AL(t)*conj(AS(t))-1i*f;

[VF,subs]=odeToVectorField(Eq1,Eq2,Eq3)

VF =

subs =

% disp(VF)

% disp(subs)

% ftotal=matlabFunction(VF);

%% Solving the ODEs

ftotal = matlabFunction(VF, 'vars', {'t', 'Y'});

ftotal = function_handle with value:

@(t,Y)[conj(Y(3)).*Y(2).*3.7292e+11i;Y(1).*Y(3).*3.7292e+11i;Y(4);conj(Y(1)).*Y(2).*3.771114367965276e+10+Y(3).*1.815225352112676e+19i+Y(4).*(-3.590000000000001e+8+1.011267605633803e+11i)+-1.035299568074342e+26-9.378225606249226e+25i]

tspan = [1 2];

ic = [0 0 r 0]; % <-- 4 states requires 4 initial values (check the order)

ySol = ode15s(ftotal, tspan, ic);

Warning: Failure at t=1.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

%% Plotting section

% plot(t, sol(:, 1), 'LineWidth', 2, 'DisplayName', 'Eq1');

% hold on;

% plot(t, sol(:, 2), 'LineWidth', 2, 'DisplayName', 'Eq2');

% plot(t, sol(:, 3), 'LineWidth', 2, 'DisplayName', 'Eq3');

% xlabel('Time');

% ylabel('Solution');

% title('Solution of the Coupled ODE System');

% grid on;