MATLAB solution to system of ODEs with forward and backward propagation

[phy127]

I have a system of coupled ODEs which tells the propagation of power P_i in an optic fiber.

<br /> \frac{\partial P_i }{\partial z} = \left (N\sigma - 1 \right ) P_i<br />
where
N = \frac{\sum_i \alpha_i P_i}{\sum_i \beta_i P_i + 1}<br />

If the signals are copropagating, there is no problem since it is easily solvable with ode45 in MATLAB. However, since there are signals which are propagating in the opposite direction, solutions must be relaxed (according to the Book) so that the forward and backward propagation powers would agree.

I'm solving the system from z = [0,L]. For those forward propagating signals, P_i(0)=finite. For those backward propagating signals, P_i(L) = some value. Can I solve this using initial value problem techniques in MATLAB? or do i need to move to boundary-value problem techniques. The thing is, I only know one side of the boundary, the value on the other side is my objective.

Please help. :)

 

[kai_sikorski]

The following simple problem has the essential features you're asking about, so let me show you how I would do it in this easier context.

\left(<br /> \begin{array}{c}<br /> x \\<br /> y \\<br /> \end{array}<br /> \right)&#039;=\left(<br /> \begin{array}{cc}<br /> -2 &amp; 1 \\<br /> 1 &amp; -2 \\<br /> \end{array}<br /> \right)\left(<br /> \begin{array}{c}<br /> x \\<br /> y \\<br /> \end{array}<br /> \right)
x(0) = 1
y(1) = 1

Here is the MATLAB code to solve it by setting up a simultaneous system for all the variables

function [ x y z] = example( N )
%Discretize the z interval with N points
z = linspace(0,1,N);
dz = 1/(N-1);
%Create the matrix system we need to solve. In the matrix the first 1 to N
%positions correspond to xs and the N+1 to 2N postions correspond to ys

A = zeros(2*N,2*N);
b = zeros(2*N,1);

%now the equations at each data point using a forward and backward
%difference for the y and x derivatives respectively
for i=1:N-1,
    %equations corresponding to x points
    A(i+1,i) = -1/dz;
    A(i+1,i+1) = 1/dz;
    A(i+1,i+1) = A(i+1,i+1)+2;
    A(i+1,i+1+N) = -1;
    %equations corresponding to y points
    A(i+N,i+N+1) = 1/dz;
    A(i+N,i+N) = -1/dz;
    A(i+N,i+N+1) = A(i+N,i+N+1) +2;
    A(i+N,i) = -1;
end
%handle end points as special cases
A(1,1) = 1; b(1) = 1;
A(2*N,2*N) = 1; b(2*N) = 1;
vars=A\b;
x = vars(1:N);
y = vars(N+1:2*N);
end

The exact solutions are
x(z) = \frac{e^{-3 z} \left(e^2-e^3+e^{2 z}+e^{3+2 z}\right)}{1+e^2}
y(z) = \frac{e^{-3 z} \left(-e^2+e^3+e^{2 z}+e^{3+2 z}\right)}{1+e^2}

If you run it you'll see it works.

Now your problem will be harder because it is non-linear which means the matrix A will actually be some Jacobian, and you'll have to do Newtons method or something like that to find the solution.

 

[kai_sikorski]

The alternative, which might be easier to code but is less elegant is to just guess the initial conditions for the variables for which you don't have them, integrate the system and then iterate until you find a set of ICs such that the right BCs are satisfied.

 

[phy127]

Thanks kai_sikorski

Your first reply is elegant! But unfortunately, my system is nonlinear. I have been searching around too and found that your second reply is the commonly used method.

On the other hand, I'm interested in the MATLAB function you presented. As for my problem, the number of simultaneous equations could be more than 50 if needed and about a quarter of that is backward propagating. I don't know yet if MATLAB can handle it in a faster rate using your second reply.

 

[kai_sikorski]

phy,

As I said in my first post you can handle non linear problems. If you're trying to solve y'=f(y) you'll need to find the Jacobian of f and do Newton's method. If there are 50 equations, and say you want to discretize using 100 z steps, that would be 5000x5000 system you'll have to solve at each iteration. It might not be fast but I think that's doable.

EDIT Sorry it won't be a jacobian of f strictly speaking but rather a function related to f. You can see the pdf file I post below.

 

[phy127]

Sorry if I misunderstood you kai. I'm not really good with linear algebra and I don't know exactly how can I do the method you are presenting. I will give it a try if you will help me lay it all. Of course that will consume your time.

 

[kai_sikorski]

This file seems to go over what you need to do

http://www.cs.rpi.edu/~flaherje/pdf/ode6.pdf

Section 6.3 first goes over a linear problem similarly to what I did, then on page 22 they talk about a non-linear problem using Newtons method.

They use centered differences instead of forward and backward like I did. You should change their expressions to use forward and backward (depending on whether you're looking at a variable that has a left or right b.c.).

Maybe there's also a way to use centered difference in your problem. I couldn't figure out how to get as many equations as unknowns that way because you run into problems when for example trying to write an equation for x_N in my post. But maybe I'm missing something and maybe someone else can chime in.

 

[kai_sikorski]

If you have Mathematica you might also be able to be lazy and just use NDSolve.

 

[phy127]

I'm reading the pdf file you directed. It's too hard for me. I don't have very good mathematical foundations. Maybe this needs a lot of time. After some time, I will post a reply if things are getting better, or worse.

About Mathematica, I'll stick with Matlab as I need to do this in Matlab or Python.

 

[kai_sikorski]

Yeah the Newton method stuff can be a little confusing.

It might be easier to try to implement a shooting method. The file I linked to has a section on that as well, or if that's too hard I'm sure you can find some other resources for the shooting method. Do a little reading and try to implement it. Once you give it a try if you have any problems go ahead and post and I think it will be easier to address specific issues than to try to give you general guidance.

Hang in there!