Analytical or numerical Solution for a coupled differential equation.

4 次查看（过去 30 天）

显示更早的评论

PS 2022-6-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1738410-analytical-or-numerical-solution-for-a-coupled-differential-equation

编辑： David Goodmanson 2022-10-10

DE.PNG

I need a help in analytical/numerical solution to a coupled differential equation as attached in the image. Where m=n=3

22 个评论
显示 20更早的评论隐藏 20更早的评论

PS 2022-10-8

编辑：Torsten 2022-10-9

在 MATLAB Online 中打开

numerical_vs_analytical_5oct.m

Thankyou @Walter Roberson and @Torsten for your replies. I have defined the span beginning and end which is [0 10], ode45 choses the step size, which is then applied in the analytical expression to compare the 2 methods. I am doing this process for couple of equations here, the coefficient named coupling_coefficient_12 is changed in the ode45 equation and the step size is chosen by the ode45 method, the doubt is that the step size used during different coupling coefficients value is significantly different, is this expected or am doing anything wrong here ? I ahve attached the code for your reference.

% Comparison between RK45 method and the analytical solution or a 2 mode

% coupled differential equation

clc;clear all;close all;

cc=0;

for Coupling_Coefficient_12=0:4:20 %coupling coefficient (not solving equation 4 and directly defining coupling coefficient values

cc=cc+1;

del_beta=0.3;

%Using numerial solution ie RK 45 to solve a coupled differential equation

%for 2 modes

options=odeset('RelTol',1e-6,'AbsTol',1e-9);

f=@(z,x) [-1i*Coupling_Coefficient_12*x(2)*exp(1i*del_beta*z);-1i*conj(Coupling_Coefficient_12)*x(1)*exp(-1i*del_beta*z)]; %x(2):LP11a, x(1):LP01

[zz,xa]=ode45(f,[0 10],[1;0],options); %[LP01=1 LP11a=0] represent the initial condition

% [0 10] is the span of z starting from 0 to 10, the steps size chosen

% by ODE45 will later be used for solving by analytical method

% xa(:,1) represent LP01 soluton

% xa(:,2) represents LP11a solution

%Analytical solution obtained using substitution elimination to solve a

%coupled differential equation or 2 mode

%Using the same "adaptive step size for z" and substituting in numerical solutions to compare the differnce between both methods

k=Coupling_Coefficient_12;

s=sqrt((Coupling_Coefficient_12*conj(Coupling_Coefficient_12))+((del_beta/2)^2));

A1_0=1; % defining initial condition LP01=1

A2_0=0; % defining initial condition LP11a=0

A1=(exp(1i*(del_beta/2)*zz')).*((cos(s*zz')-(1i*(del_beta/2)*sin(s*zz'))/s).*A1_0-(((1i*Coupling_Coefficient_12*sin(s*zz')/s)).*A2_0));

A2=(exp(-1i*(del_beta/2)*zz')).*((-1i*conj(Coupling_Coefficient_12)*(sin(s*zz'))/s).*A1_0+(cos(s*zz')+1i*(del_beta/2)*(sin(s*zz'))/s).*A2_0);

diff_method1(cc)=mean(abs(xa(:,1)-A1.').^2); %difference between LP01 solution obtained numerically vs analytically

diff_method2(cc)=mean(abs(xa(:,2)-A2.').^2); %difference between LP11a solution obtained numerically vs analytically

figure(cc)

plot(zz,xa)

end

Warning: Imaginary parts of complex X and/or Y arguments ignored.

figure(cc+1);

plot(0:4:20,diff_method1)

hold on

plot(0:4:20,diff_method2)

Torsten 2022-10-9

编辑：Torsten 2022-10-9

Analytical solution and solution from ODE45 are almost identical. So it seems ODE45 needs these smaller time steps to get the correct solution within the prescribed error tolerance.

The coupling coefficient is directly related to the frequency of the sin and cos terms in the analytical solution. And if you plot sin(x) and sin(20*x), you will see that it will be much more difficult to resolve the cycles of the trigonometric functions for bigger coupling coefficients.

Plotting the solutions A1 and A2 might help in understanding why solutions for bigger coupling coefficients are more complicated to get (see above).

David Goodmanson 2022-10-10

编辑：David Goodmanson 2022-10-10

在 MATLAB Online 中打开

Hi PS,

If the elements of Cmn are differing functions of z, then an analytic solution will be rare. If all the Cmn are constants, the following is a solution where C is a square matrix of arbitrary size, within reason.

Let B be the diagonal matrix whose k,k element is beta_0k. Then for the matrix exponential

expB(z) = expm(i*B*z)     % diagonal matrix
exp(i*beta_0k*z)          % its k,k th element

In matrix notation the original equation is

dA/dt = -i*expB(z)*C*expB(-z)*A,

where A is a column vector with n components.

For the solution, let

[V lambda] = eig(B+C)

Lambda is the diagonal matrix of eigenvalues, and denoting its matrix exponential in a similar way as before,

expL(-z) = expm(-i*lambda*z).         % diagonal matrix
exp(-i*lambda_k*z)                    % its k,k the element

The solution is

A = expB(z)*V*expL(-z)*g

where g is a column vector of constant amplitudes that are determined by initial conditions. If the initial conditions are set as A = A0 at z = 0 for some column vector A0, then

g = V\A0.

请先登录，再进行评论。

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