Including a periodic piecewise function of time in coupled ODE

Hi All,

I want to solve an ODE which includes a function which is periodic and peicewise contructed (although not with the peicewise functionality...). This is something very similar to:

https://uk.mathworks.com/matlabcentral/answers/97074-how-do-i-solve-an-ode-with-time-dependent-parameters-in-matlab

However not the same, I think.

I have made a go at trying to have the function defined symbolically with in its period and as just a general function for a number of periods. I then show how I solve the coupled ODE's without this function as a coefficient.

I have then commented out the code which is my (wrong) way of trying to incorporate the function into the ODE45 so that the code will run.

At the end I have plotted my solution for the working ODE without the function and the function itself to demonstrate that both are solved over the same x-domain.

The peicewise function is in black, and as you can see it is smooth, continuoud and differntiable.

I am not adverse to redefining the function if it can be done, or smoothing or interpolating the function in the ODE solver.

clear
a=0.15;
b=0.2;
m=448;


x=linspace(0,m,10000);

intvl = [0 3*m];

eta= @(x) [(0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))))];



etafull = repmat(eta(x),1,3);
xfull=linspace(intvl(1), intvl(2),length(etafull));

figure(1)
plot(xfull,etafull)
xlim([0 2.5*m])


w0 = 2*pi*67*10^9;

wj = 2*pi*86*10^9;

wp = 2*pi*12*10^9;

wsfac = 0.6;
wifac = 1-wsfac;
        
ws = wsfac.*wp;
w2p = 2*wp;
Ap0 = 0.5*w0/wp;
As0 = Ap0*sqrt(0.0057*wp/ws);
%As0
A2p0 = 0;
wi = wifac.*wp;
kp = (wp/w0)*(1/(sqrt(1-(wp/wj)^2)));
ks = (ws/w0)*(1/(sqrt(1-(ws/wj)^2)));
ki = (wi/w0)*(1/(sqrt(1-(wi/wj)^2)));
k2p = (w2p/w0)*(1/(sqrt(1-(w2p/wj)^2)));

delk = 3*ws*wi*wp/(2*w0*(wj^2));

modk = sqrt(wp.*Ap0^2/(ws*As0^2 + wp.*Ap0^2));

maxBeta=0.433

dA = @(xfull,A)[-(maxBeta/2)*ks*ki*A(2)*A(3)*exp(1i*(ks+ki-kp)*xfull) + (maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*xfull);
(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*(kp-ki-ks)*xfull);
 (maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*(kp-ks-ki)*xfull);
 -(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*xfull)];

[xfull,A] = ode45(dA, xfull ,[Ap0; As0; 0; 0]);

% dA = @(x,A)[-eta*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
% eta*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
%  eta*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
%  -eta*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];
% 
% [x,A] = ode45(dA, x ,[Ap0; As0; 0; 0]);


P=[(wp.^2)*(abs(A(:,1))).^2/((wp.^2)*(abs(A(1,1))).^2), (ws.^2)*(abs(A(:,2)).^2)/((wp.^2)*(abs(A(1,1))).^2), (wi.^2)*(abs(A(:,3)).^2)/((wp.^2)*A(1,1).^2), (w2p.^2)*(abs(A(:,4)).^2)/((wp.^2)*A(1,1).^2)];

figure(2)
plot(xfull,P)
hold on
plot(xfull,etafull,'k','LineWidth',3)
hold off


% dA = @(x,A,eta)[-eta*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
% eta*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
%  eta*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
%  -eta*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];
% 
% [x,A] = ode45(dA, x ,[Ap0; As0; 0; 0]);

Post edit: the variable x is the same variable to solve for A and that eta is defined over, that is and are over the same domain/variable x.

Thank you for any help,

Tom

Answer 1

[Alan Stevens]

Because eta is itself a function of x, you need to have

 dA = @(x,A)[-eta(x)*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
 eta(x)*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
  eta(x)*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
  -eta(x)*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];

[Thomas Dixon]

Thanks a lot Alan.

[Thomas Dixon]

Hi Alan,

Suppose I wanted to extend the range of this solution to the range of 'etafull' (or more), How would I go about this. 'eta' is defined over the range of the period x = [0 m] but actually I may want to do this for a range [0 3*m]. Where I do this above I use the 'repmat' function to repeat the function but now it loses its symbolic reference to x, ie it wont work in the ODE solver.

Is there a way to repeat the function periodically and have matlab give me the x dependence (essentially this just means replacing the x --> x-N*T at each repetition of the function) where N is the number of periods and T is the period.

If not is there a way to sample the function at the apropriate places for the ODE solver, or something similar.

In the end I would like to create some function which is peicewise defined. I would like to be able to repeat this function (that is repeat the y-values periodically, not repeat the function y(x) such that as x grows with each period y(x+T) =/= y(x)) such that it can be included in an ODE solver. I would then like to be able to solve that ODE across any domain size (I accept differentiability, tolerances etc may be a problem but I mean that so long as the peicewise function I have defined is acceptable it has a go rather than just returning an error about arguments etc.).

I hope this makes sense,

Tom

[Alan Stevens]

What you want should work ok if you have:

dA = @(x,A)[-eta(x).*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) ...
    + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
    eta(x).*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
    eta(x).*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
    -eta(x).*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];

[x,A] = ode45(dA, [0 3*m] ,[Ap0; As0; 0; 0]);

However, ode45 can't meet the integration tolerances at around x = 353. This suggests you should take a very careful look at your (complicated) function, and/or your data values.

[Thomas Dixon]

Thanks again Alan,

I have implemented this. But I find that my function is not defined for these limits. This code below plots the function eta(black) and the solutions on top of each other and because the function is now not defined It stops after on cycle of itself.

clear
a=0.15;
b=0.2;
m=448;


x=linspace(0,m,10000);

intvl = [0 3*m];

eta= @(x) [(0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))))];



etafull = repmat(sym(eta(x)),1,3);
xfull=linspace(intvl(1), intvl(2),length(etafull));

figure(3)
plot(xfull,etafull)
xlim([0 2.5*m])



Cj=235e-15;
Ic=1.8e-6;
maxBeta=0.25;
C0=32e-15;
L=88e-12;
w0=1/sqrt(L*C0)
wj=1/sqrt(L*Cj)
wp = 2*pi*15*10^9;

Ap0=0.53*w0/wp;

% w0 = 2*pi*67*10^9;
% 
% wj = 2*pi*86*10^9;


wsfac = 0.6;
wifac = 1-wsfac;
        
ws = wsfac.*wp;
w2p = 2*wp;
% Ap0 = 0.5*w0/wp;
As0 = Ap0*sqrt(0.0057*wp/ws);
%As0
A2p0 = 0;
wi = wifac.*wp;
kp = (wp/w0)*(1/(sqrt(1-(wp/wj)^2)));
ks = (ws/w0)*(1/(sqrt(1-(ws/wj)^2)));
ki = (wi/w0)*(1/(sqrt(1-(wi/wj)^2)));
k2p = (w2p/w0)*(1/(sqrt(1-(w2p/wj)^2)));

delk = 3*ws*wi*wp/(2*w0*(wj^2));

modk = sqrt(wp.*Ap0^2/(ws*As0^2 + wp.*Ap0^2));

maxBeta=0.433

dA = @(x,A1)[-(maxBeta/2)*ks*ki*A1(2)*A1(3)*exp(-1i*(delk)*x) + (maxBeta/2)*k2p*A1(4)*kp*conj(A1(1))*exp(1i*(k2p-2*kp)*x);
(maxBeta/2)*ki*kp*conj(A1(3))*A1(1)*exp(1i*(delk)*x);
 (maxBeta/2)*ks*kp*conj(A1(2))*A1(1)*exp(1i*(delk)*x);
 -(maxBeta/4)*kp^2*A1(1)^2*exp(1i*(2*kp-k2p)*x)];

[x1,A1] = ode45(dA, [0 m*3] ,[Ap0; As0; 0; 0]);

P1=[(wp.^2)*(abs(A1(:,1))).^2/((wp.^2)*(abs(A1(1,1))).^2), (ws.^2)*(abs(A1(:,2)).^2)/((wp.^2)*(abs(A1(1,1))).^2), (wi.^2)*(abs(A1(:,3)).^2)/((wp.^2)*A1(1,1).^2), (w2p.^2)*(abs(A1(:,4)).^2)/((wp.^2)*A1(1,1).^2)];


% dA = @(xfull,A)[-etafull(xfull)*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*xfull) + etafull(xfull)*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*xfull);
% etafull(xfull)*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*xfull);
%  etafull(xfull)*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*xfull);
%  -etafull(xfull)*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*xfull)];
% 
% [x,A] = ode45(dA, xfull ,[Ap0; As0; 0; 0]);




dA = @(x,A2)[-eta(x)*(maxBeta/2)*ks*ki*A2(2)*A2(3)*exp(-1i*delk*x) + eta(x)*(maxBeta/2)*k2p*A2(4)*kp*conj(A2(1))*exp(1i*(k2p-2*kp)*x);
eta(x)*(maxBeta/2)*ki*kp*conj(A2(3))*A2(1)*exp(1i*delk*x);
 eta(x)*(maxBeta/2)*ks*kp*conj(A2(2))*A2(1)*exp(1i*delk*x);
 -eta(x)*(maxBeta/4)*kp^2*A2(1)^2*exp(1i*(2*kp-k2p)*x)];

[x2,A2] = ode45(dA, [0 m*3] ,[Ap0; As0; 0; 0]);



P2=[(wp.^2)*(abs(A2(:,1))).^2/((wp.^2)*(abs(A2(1,1))).^2), (ws.^2)*(abs(A2(:,2)).^2)/((wp.^2)*(abs(A2(1,1))).^2), (wi.^2)*(abs(A2(:,3)).^2)/((wp.^2)*A2(1,1).^2), (w2p.^2)*(abs(A2(:,4)).^2)/((wp.^2)*A2(1,1).^2)];

figure(4)
plot(x1,A1,'b')
hold on
plot(x2,A2)
fplot(eta,'k','LineWidth',3)
xlim([min(x1) max(x1)])
ylim([-1.5 10])
hold off

figure(5)
plot(x1,P1,'b')
hold on
plot(x2,P2)
fplot(eta,'k','LineWidth',3)
xlim([min(x1) max(x1)])
ylim([-1.1 1.1])
hold off

Note there are some small changes to my initial conditions which makes each solution fail at separtate points to yours but when I plot(x1,A1) it gives me the entire 3*m solution even though it only reports giving me it up to x=765.

This has thrown me a little bit

[Thomas Dixon]

Ah the x vector is spaced non-linearly. I have not seen this before. Is this something unphysical? or is it trying to give me something based on the tolerances. I think I can see that becuase eta has gone to zero as it is undefined then it is obvious intuitively that the derivative is now zero everywhere and thus it doesnt change. But is MATLAB clever enough to know this by 'looking' at it symbolicaly or is it solving it regulalry and saying 'ah well this solution is within a change of y that is aceptable so dont bother solving it at smaller x-step'. Also this behaviour is obviously not what I want.

I would like to be able to periodicaly repeat the eta function over the x-domain of the ODE so that the eta(x) function does not go to zero. Then solve the ODE. I don't mind defining this domain outside the ODE function previosuly but I would like it if I didnt have to manually write down the new function for each period of oscillation.

[Walter Roberson]

piecewise functions can only be used with ode*() functions if they have C2 or higher continuity (their second derivatives must be continuous.) The ode*() routines will typically not detect the problem if the functions are C1 continuous (first derivative is continuous) and will instead silently produce the wrong value. The ode*() routines will typically not detect problems with C1 continuity if the jump is proportionately less than the current integration tolerance, and again will silently return the wrong value in such a case.

[Walter Roberson]

Ah the x vector is spaced non-linearly. I have not seen this before. Is this something unphysical?

The ode*() functions do not return values at constant time (or in your case, x) intervals unless you pass in a tspan that has more than 2 elements and the tspan is spaced regularly.

The ode*() functions are all adaptive functions, changing the step interval as needed to meet local conditions. In places where the integration is going well, it will increase the step size. It will keep increasing the step size as long as it is successful in doing so. The functions make calculations up to a certain order, and use the information to deduce a new state. They then use that deduced state to predict at one or more additional locations, and evaluate at those locations and cross-check the actual values against predictions; if the result is within tolerance then the deduced state is "accepted" (and next round, the step size will be increased.) If the actual values were not within tolerance of those predicted by the deduced state, then the trial step is rejected, and the step size is decreased and a new series of measurements is taken relative to the previous state. So as long as everything goes well, the step size accelerates, and when it figures out that it has overstepped, it retreats back to a smaller step.

This is why the x vector is spaced non-linearly: the spacing used is adaptive, more sparse when that works, less sparse when the conditions are changing rapidly.

[Thomas Dixon]

Thanks Walter. I will look to differentiate the function and check it has a continuouse second derivative.

Did you try plotting the function I defined? As you can see it oscilates between -1 and 1 on a period of x=433.

Is it possible to repeat this function periodically and define at as a new function but transfer over its x dependence? such that this new function could be put into an ode.

Thats is turn eta(x) into some function eta'(x) where eta'(x) is say 433*3 = 1299 long and then MATLAB knows that each position is for x=1,2,3.... 1299. Then when it goes into the ODE it can linearly interpolate between points to get an approximate value etc?

As it is a coefficient I am hoping to solve the ODE for A(x) but I just need to assign a value to eta(x) at each point its exact value is not so important right now, hopefully.

[Walter Roberson]

In particular you would need to match the derivatives at the boundaries.

If the derivatives at the boundaries are matched, in some cases you can use argument reduction, if the periods are all the same length and are identical. If, though, the duration stretches as you go, or the amplitude changes as you go, argument reduction alone might not work.

Linear interpolation is not C1 continuous; you need spline for C1 and C2 continuity.

[Thomas Dixon]

When you say the boundaries do you mean the bondaries of the integration or the boundaries of each section of peicewise function?

The peicewise function eta is meant to be continuous as you can see from the plot and so shouldn't be discontinuous at the boundaries.

[Walter Roberson]

boundaries of each piecewise section.

[Thomas Dixon]

I am very grateful for the depth of this, I think I may have confused the problem.

I want to solve some coupled differential equation of the form:

dA(1)/dx = -eta[x]*A(2)[x]

dA(2)/dx = eta[x]*A(1)[x]

The integration function I want to find is A, and while eta[x] is a function of x I wouldn't care if the value given to the ODE solver is just a value of eta defined at certain x points (these points are then defined by where the ODE solver is numerically computing the derivative), then the ODE solver treats eta as just a number which it may/may not need to find depending on where it wants to solve the ODE. (I hope that makes sense, I do understand that a constant that changes is not a constant, but I am looking for insight into how this may be done without defining it in some crude made up way from my head)

But to be clear the the function eta is meant to be continuous, differentiable and single valued, as shown in the plot:

I believe it is very similar to a two sided bump function which has no analytic form for its entire period (I am happy to be proven wrong). I may then want to extend this (see eta full) which would then definitely not have an analytic form (Very happy to be proven wrong).

The ODE solver works very well in the period of eta(x) this is because eta(x) is defined in that period. However then eta(x) goes to 0 because it hasnt been repeated (that is eta_full). but now eta_full is not a function of x and so does not go into the ODE solver. The reason the ODE solver can comlete solving though is because MATLAB has defined the peicewise function eta to go to zero outside of its period. This makes solving the differential equation very easy as now al derivatives are proportional to zero, the function doesnt change and you can see the reuslt in this plot:

Where the blue circles are solutions when eta is not included and the coloured lines are when eta is included. Note the integrator did not 'fail' it was the combination of an ill defined function as a coefficient, an illdefined ODE and my superior lack of pragramming ability, I believe this can be sorted but I am not making my intentions clear enough.

I appreciate your continued support,

Tom

Ps this is obviously going to be extended so I would like a solution where I don't have to define the peicewise function for the entire duration by stitching it together multuiple times so some repeated functionality would be desirable.

[Walter Roberson]

Argument reduction.

eta(ReduceEtaArgument(x))

where ReduceEtaArgument(x) is something like

mod(x, duration_of_one_eta_cycle)

[Alan Stevens]

You can repeat your eta function shape for as long as you like using the following function

function etaval = eta(x)

        a=0.15;
        b=0.2;
        m = 448;
        x = mod(x,m);
        
        etaval=  (0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m)) ...
            +exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) ...
        + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))...
        +exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))));

end

Your dA function still gives problems around x = 352 though.

[Thomas Dixon]

Thanks both.

Alan - Do you reproduce the figure I put in my comment above. Here the dA function has worked beyond 352 and for its entire range. I would like to know why mine has worked and yours hasn't. BE AWARE IT IS HIGHLY LIKELY THAT I EDITED MY DA FUNCTION IN ORIGINAL POST DURING OUR DISCUSSION AND SO YOU MAY HAVE A DIFFERENT ONE TO THE ONE I AM USING.

Walter - I think your solution is essentially the same as Alan's which I have somehow got working. If you disagree then let me know why, this is very interesting to me. Thank you for such detailed response.

Both - Somewhere on my PC in the same fodler as this script i have saved the function:

function etaval = dAeta(x)
        a=0.15;
        b=0.2;
        m = 448;
        x = mod(x,m);
        
        etaval=  (0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m)) ...
            +exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) ...
        + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))...
        +exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))));
end

With the name dAeta.m

I then run this code:

clear
a=0.15;
b=0.2;
m=448;

x=linspace(0,3*m,10000)
intvl = [0 3*m];

%eta= @(x) [(0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))))];

x=mod(x,m)

etafull = repmat(sym(eta(x)),1,3);
xfull=linspace(intvl(1), intvl(2),length(etafull));


Cj=235e-15;
Ic=1.8e-6;
maxBeta=0.25;
C0=32e-15;
L=88e-12;
w0=1/sqrt(L*C0)
wj=1/sqrt(L*Cj)
wp = 2*pi*15*10^9;

Ap0=0.53*w0/wp;


wsfac = 0.6;
wifac = 1-wsfac;
        
ws = wsfac.*wp;
w2p = 2*wp;
As0 = Ap0*sqrt(0.0057*wp/ws);
A2p0 = 0;
wi = wifac.*wp;
kp = (wp/w0)*(1/(sqrt(1-(wp/wj)^2)));
ks = (ws/w0)*(1/(sqrt(1-(ws/wj)^2)));
ki = (wi/w0)*(1/(sqrt(1-(wi/wj)^2)));
k2p = (w2p/w0)*(1/(sqrt(1-(w2p/wj)^2)));

delk = 3*ws*wi*wp/(2*w0*(wj^2));

modk = sqrt(wp.*Ap0^2/(ws*As0^2 + wp.*Ap0^2));

maxBeta=0.433

dA = @(x,A1)[-(maxBeta/2)*ks*ki*A1(2)*A1(3)*exp(1i*(ks+ki-kp)*x) + (maxBeta/2)*k2p*A1(4)*kp*conj(A1(1))*exp(1i*(k2p-2*kp)*x);
(maxBeta/2)*ki*kp*conj(A1(3))*A1(1)*exp(1i*(kp-ki-ks)*x);
 (maxBeta/2)*ks*kp*conj(A1(2))*A1(1)*exp(1i*(kp-ks-ki)*x);
 -(maxBeta/4)*kp^2*A1(1)^2*exp(1i*(2*kp-k2p)*x)];

[x1,A1] = ode45(dA, [0 m*3] ,[Ap0; As0; 0; 0]);

P1=[(wp.^2)*(abs(A1(:,1))).^2/((wp.^2)*(abs(A1(1,1))).^2), (ws.^2)*(abs(A1(:,2)).^2)/((wp.^2)*(abs(A1(1,1))).^2), (wi.^2)*(abs(A1(:,3)).^2)/((wp.^2)*A1(1,1).^2), (w2p.^2)*(abs(A1(:,4)).^2)/((wp.^2)*A1(1,1).^2)];

dA = @(x,A2)[-eta(x)*(maxBeta/2)*ks*ki*A2(2)*A2(3)*exp(1i*(ks+ki-kp)*x) + eta(x)*(maxBeta/2)*k2p*A2(4)*kp*conj(A2(1))*exp(1i*(k2p-2*kp)*x);
eta(x)*(maxBeta/2)*ki*kp*conj(A2(3))*A2(1)*exp(1i*(kp-ki-ks)*x);
 eta(x)*(maxBeta/2)*ks*kp*conj(A2(2))*A2(1)*exp(1i*(kp-ks-ki)*x);
 -eta(x)*(maxBeta/4)*kp^2*A2(1)^2*exp(1i*(2*kp-k2p)*x)];

[x2,A2] = ode45(dA, [0 3*m] ,[Ap0; As0; 0; 0]);

P2=[(wp.^2)*(abs(A2(:,1))).^2/((wp.^2)*(abs(A2(1,1))).^2), (ws.^2)*(abs(A2(:,2)).^2)/((wp.^2)*(abs(A2(1,1))).^2), (wi.^2)*(abs(A2(:,3)).^2)/((wp.^2)*A2(1,1).^2), (w2p.^2)*(abs(A2(:,4)).^2)/((wp.^2)*A2(1,1).^2)];

figure(4)
plot(x1,A1,'bo')
hold on
plot(x2,A2)
%fplot(eta,'k','LineWidth',3)
xlim([min(x1) max(x1)])
ylim([-1.5 10])
hold off

figure(5)
plot(x1,P1,'bo')
hold on
plot(x2,P2)
% fplot(eta,'k','LineWidth',3)
%fplot(x,eta)
xlim([min(x1) max(x1)])
ylim([-1.1 1.1])
hold off


%%%%%%---------------------------------------------%%%%%%%%
% 
% function etaval = dAeta(x)
%         a=0.15;
%         b=0.2;
%         m = 448;
%         x = mod(x,m);
%         
%         etaval=  (0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m)) ...
%             +exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) ...
%         + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))...
%         +exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))));
% end
% 

Which somehow works, and gives me a sensible plot. The sensibility is quite nice actually if you notice the turning points of figure 4 align with the half period of the eta function such that the derivatives 'no longer turn' when you include the eta function, this is the behaviour i want.

I now have a much nicer set of problems:

- I dont understand how it knows what eta is

- I dont know how to plot eta, before I used fplot(eta) but now that says I dont have the input arguments

This is very exciting times, thank you

Tom

[Thomas Dixon]

Aplogies here are the plots

[Thomas Dixon]

Ignore that last bit it is a fundemental misunderstanding of what a function and input arguments are.

I guess the ODE solver also goes to this function at any x point it wants and asks 'what the value is' arbitrarily.

It is then simple I assume to pass a vectorised x (ie the one each ode solver returns me x1, x2) to the same function and simply plot the result that is:

eta(0,m/2,m) = [0,0,0]

eta(0:1:m/2) = ['whatever numbers make that tapered sine wave evaluated at x=0:1:m/2']

so I can always ask for

eta(some_vector)

and then:

plot(some_vector,eta)

to see eta.

Lovely