Analytical or numerical Solution for a coupled differential equation.

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

[Torsten]

What is t ? A second independent variable ?

[PS]

Thankyou for your reply. This is a coupled mode differential equation for the modes propagating in a perturbed fiber. For m=n=3, ignoring t (time) for now. I have elaborated the equations. There are 3 equations which needs to be solved (attached Equation.png). Cmn is a 3x3 matrix dependent on z and beta_m - beta_n are constant values. Z is varying in small steps and i need the solution of A1,A2,A3 in terms of C,beta,z. For 2 mode equation i could solve it analytically but for 3 modes the solution is becoming tedious. I need help in solving these 3 equations either analytically or numerically. Since we have complex i term, so i am not sure which method could be used or if RK method can be helpul here.

Thanks in advance

[Torsten]

Should be easy to solve with ODE45. The solver can handle complex-valued ordinary differential equations.

Give it a try.

[PS]

Any idea if there exists an analytical solution for this ?

[Torsten]

The only thing you could try is MATLAB's "dsolve" or MATHEMATICA for an analytical solution.

[PS]

Hey, I tried solving using dsolve, but i get a "Warning: Unable to find explicit solution".

[Torsten]

Then - as said - use a numerical ode solver, e.g. ode45.

[David Goodmanson]

Hi PS,

For the two mode case was there a specific form of Cmn(z) that allowed an analytic solution? That won't be the case in general.

[PS]

Hi David,

Thanks for your reply, Cmn(z) is a variable dependent on z, which comes after solving an integral. While solving the dsolve i used it as a constant (eg Cmn=2) and tried solving the equation but I still got the warning. I guess dsolve doesnt solve for complex equations,ryt ?

[Walter Roberson]

dsolve can handle complex-valued equations. However, it is not all that good with complicated equations.

A trick with dsolve is that if it cannot handle the equations with initial conditions, then sometimes it can handle them without initial conditions, and then you might be able to guide it through to solve for the generated constants. That said, there are a lot of cases that dsolve cannot handle even without initial conditions.

[PS]

Hi,

I tried solving the 2 mode equation using RK45 in matlab. I have defined the relative error tolerance of 10^(-6) and absolute error tolerance of 10^(-9) in matlab while using ODE45 function, and the step size to be 10^5. Can we define both error tolerance and step size together? If not, how do we find the step size for the defined error tolerance ?

[Torsten]

ODE45 calculates its own step size - you can't prescribe it.

You can only choose a maximum and a minimum stepsize.

You can define all quantities in the options structure independently from each other.

[PS]

Thankyou for your reply, actually i was trying to compare 2 mode solution using ODE45 with analytical solution available with me, to understand the error between the two. As you mentioned ODE45 calculates its own step size, is there a way to know that ? As i will be giving the same step size in my analytical method to have a fair comparison between the 2 techniques.

[Torsten]

As you mentioned ODE45 calculates its own step size, is there a way to know that ?

You can get an impression of the stepsize from ODE45 by choosing

tspan = [tstart tend]

as a two-element vector.

Then the output T and Y will be given for the (successful) integration steps of ODE45:

[T, Y] = ode45(...) 

As i will be giving the same step size in my analytical method to have a fair comparison between the 2 techniques.

An analytical method does not need a stepsize.

[Walter Roberson]

The ode* functions are all adaptive step size. Each time a step succeeds they increase the step size, and each time a step fails they decrease it. If you set the option Refine to 1 and your tspan is exactly two elements then an output will be recorded for each step that succeeds, so you could diff(t) to see the instantaneous time step.

[Torsten]

@PS

Walter is right - additionally to the two-element vector for tspan, you have to set Refine to 1 when using ODE45.

This will give you the solution at the end of all successful time steps.

[Walter Roberson]

And of course you can take the t values output by ode45 and subs() them into the dsolve solution to get the analytic solution at the same times, suitable for plotting together with a dashed or dotted line for the second line (to be able to see underneath it in places of overlap)

[PS]

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.

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

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

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

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]

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 ?

You mean that ode45 takes different step sizes for different values of "Coupling_Coefficient_12" ? This might be the case. At least the errors in both components show that the ode45 solution approximates the analytical solutions quite well, doesn't it ?

[PS]

Hi @Torsten, thankyou for your reply. Yes, ode45 is taking different step size for different coupling coeffificients which may be fine, but i was wondering if the change in step size between the different values would change so significantly as i see in my case.

Couping_coefficient_12 = 0, step size = 41

Couping_coefficient_12 = 4, step size = 929

Couping_coefficient_12 = 8, step size = 1849

Couping_coefficient_12 = 12, step size = 2773

Couping_coefficient_12 = 16 step size = 3689

Couping_coefficient_12 = 20, step size = 4609

[Torsten]

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]

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.