matlab solving numerical differential equations [pic]

192 views Asked by At

I'm trying to solve this numerical differential equations, can someone help? pic of the differential equations

clc;
clear;
syms A1(z) A2(z)
lamda1 = 1560*(10^-9);
c=3*(10^8);
d_eff=27*(10^-12);
omga1=(2*pi*c)/(lamda1);
omga2=omga1*2;
n=2.2;
k1=(n*omga1)/c;
k2=(n*omga2)/c;
ode1 = diff(A1) == (2*i*(omga1^2)*d_eff*A2*conj(A1)*exp(-i*(2*k1-k2)*z))/(k1*(c^2));
ode2 = diff(A2) == (i*(omga2^2)*d_eff.*(A1.^2).*exp(i*(2*k1-k2)*z))/(k2*(c^2));
odes = [ode1; ode2];
cond1 = A1(0) == 1;
cond2 = A2(0) == 0;
conds = [cond1; cond2];
M = matlabFunction(odes)
sol = ode45(M(1),[0 20],[2 0]);
1

There are 1 answers

7
John Bofarull Guix On

in this question both ODE are coupled, hence there's only 1 ODE to solve:

1.- use 1st equation to write

A1=f(A2,dA2/dz)

and feed this expression into 2nd equation.

2.- regroup

n1=1j*k1^4/k2^2*1/deff*.25

n2=1j*3*(k1-k2)

now the ODE to solve is

y'=n1/y^2*exp(n2*z)

3.- It can obviously be done in MATLAB, but for this particular ODE in my opinion the Wolfram online ODE Symbolic Solver does a better job.

Input the obtained ODE of previous point into the ODE solver available in this link

https://www.wolframalpha.com/input?i=y%27%27+%2B+y+%3D+0

and solve

4.- The general (symbolic) solutions for A2 are

enter image description here

Note that I used k1 instead of n1 and k2 instead of n2 just in the Wolfram ODE Solver.

Rewording ; the k1 k2 expressions of the general solutions are not the wave numbers k1 k2 of the equations in the question. Just replace accordingly.

5.- Now get A1 using expression in point 2.

6.- I have spotted 2 possible errors in the MATLAB code posted in the question that shouldn't be ignored by the question originator:

In the far right side of the posted MATLAB code, the exponential expressions

6.1.- both show

exp(-i(2*k1-k2)z))/(k1(c^2))

there's this (k1*(c^2)) dividing the exponent wheras in the question none of the exponentials show such denominator their respective exponents.

6.2.- the dk or delta k expression in the exponentials of the question are obviously k2-k1 or k1-k2 , here there may be room for a sign ambiguity, that may shoot a wave solution onto the opposite direction, yet the point here is where

*exp(-1i*(2*k1-k2)*z)

should probably be

*exp(-1i*2*(k1-k2)*z)

or just

exp(-1i*(k1-k2)*z)

6.3.- and yes, in MATLAB (-1)^.5 can be either expressed with 1j or 1i but as written in the MATLAB code made available in the question, since only a chunk of code has been made available, it's fair to assume that no such i=1j has been done.