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 .

1 次查看（过去 30 天）
Yogesh
Yogesh 2024-1-3
评论： Yogesh 2024-1-3
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.

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

采纳的回答

Sam Chak
Sam Chak 2024-1-3
Ran in:
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;

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by