Parameter variation in BVP shooting method

Hi everybody!

I wrote a code that solves the following problem:

This is a 2nd order ode and the dependant variable here is rho (tilde) and the independant variable is x (tilde). K (tilde) and B (tilde) are both parameters. The boundary conditions are:

The code is:

clear all
clc
%Shooting Method for Non-Linear ODE
k=3.25e8; %Debye length - needs to be calculated.
D=4e-9; %separation distance
y0=1; %i.c
ICguess_on_target=fzero(@(x) bar_res(x), 0.25);
[x,y]=ode45(@bar_density,[0,k.*D],[y0 ICguess_on_target]);
plot(x,y(:,1));
xlabel('\kappaD');
ylabel('$\tilde{\rho}$','Interpreter','latex')
title('$\tilde{\kappa}=\tilde{B}=1$','Interpreter','latex'); 

function dydx=bar_density(x,y)
kdim=1;Bdim=1;
dydx=[y(2);2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim];
end

function r=bar_res(IC_guess)
y0=1;
yf=1;
k=3.25e8;
D=4e-9;
[x,y]=ode45(@bar_density,[0,k.*D],[y0 IC_guess]);
r=y(end,1)-yf
end

This works just fine and I'm getting great results and plots of rho(tilde) as a function of x(tilde). Now I want to run the code in a loop and vary the parameters K (tilde) and B (tilde). For instance I would like to get 6 plots of rho(tilde) as a function of x(tilde) for 6 different sets of K (tilde) and B (tilde). I have been trying to insert a for loop inside the code but I'm obviously doing it wrong.

I would appreciate some help guys.

Thanks,

Roi

Answer 1

[Alan Stevens]

The following works:

% Shooting Method for Non-Linear ODE
k=3.25e8; %Debye length - needs to be calculated.
D=4e-9; %separation distance
y0=1; %i.c

kdimarray = 1:0.2:2; Bdimarray = 1:0.5:3.5;
ICguess_on_target = zeros(size(kdimarray));
for i = 1:6
    kdim = kdimarray(i); Bdim = Bdimarray(i);
    ICguess_on_target(i) = fzero(@(x) bar_res(x,kdim,Bdim), 0.25);
    [x,y]=ode45(@bar_density,[0,k.*D],[y0 ICguess_on_target(i)],[], kdim, Bdim);
    plot(x,y(:,1));
    hold on
end
xlabel('\kappaD');
ylabel('$\tilde{\rho}$','Interpreter','latex')
title('$\tilde{\kappa}=\tilde{B}=1$','Interpreter','latex'); 

function dydx=bar_density(~, y, kdim, Bdim)

dydx=[y(2);2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim];
end

function r=bar_res(IC_guess,kdim,Bdim)
y0=1;
yf=1;
k=3.25e8;
D=4e-9;
[~,y]=ode45(@bar_density,[0,k.*D],[y0 IC_guess], [], kdim, Bdim);
r=y(end,1)-yf;
end

[Roi Bar-On]

Thank you!! It is working great!

By the way, how can I attribute corresponding legend to each line? And eventually save the rho(tilde) values for each set (get their output)?

[Alan Stevens]

The following does it (it seems a little clumsy and there might well be a better way!)

% Shooting Method for Non-Linear ODE
k=3.25e8; %Debye length - needs to be calculated.
D=4e-9; %separation distance
y0=1; %i.c

kdimarray = 1:0.2:2; Bdimarray = 1:0.5:3.5;
ICguess_on_target = zeros(size(kdimarray));

for i = 1:6
    kdim = kdimarray(i); Bdim = Bdimarray(i);
    ICguess_on_target(i) = fzero(@(x) bar_res(x,kdim,Bdim), 0.25);
    [x,y]=ode45(@bar_density,[0,k.*D],[y0 ICguess_on_target(i)],[], kdim, Bdim);
    plot(x,y(:,1));
    
    ylen = length(y(:,1));
    Y(1:ylen,i) = y(:,1);  % saves y(:,1) values 
    hold on
    lg = ['kdim = ',num2str(kdim),' Bdim = ', num2str(Bdim)]; % prepares legend
    lgnd(i,1:length(lg)) = lg;
end
xlabel('\kappaD');
ylabel('$\tilde{\rho}$','Interpreter','latex')
title('$\tilde{\kappa}=\tilde{B}=1$','Interpreter','latex'); 
legend(lgnd)

function dydx=bar_density(~, y, kdim, Bdim)

dydx=[y(2);2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim];
end

function r=bar_res(IC_guess,kdim,Bdim)
y0=1;
yf=1;
k=3.25e8;
D=4e-9;
[~,y]=ode45(@bar_density,[0,k.*D],[y0 IC_guess], [], kdim, Bdim);
r=y(end,1)-yf;
end

[Roi Bar-On]

Hey Alan

First of all, thanks!

I think that the legend that came out is actually in the opposite order (according to my physical understanding of the problem. e.g. I know that for K=B=1 I should get the horizontal line). It seems like the legend order has flipped:(

Am I right? I'm asking it purely by technical/syntax wise.

[Alan Stevens]

The legend is in the right order. When I ran K = B = 1 on its own I didn't get a horizontal line, I got a parabolic type shape!

[Roi Bar-On]

You are right. What happend was that I wrote the euqation incorrectly. It shuold be:

Thank you!!

Roi

[Alan Stevens]

I see. Also, given the shape of the curves it might be better to put a negative initial guess for the slope in fzero.

[Roi Bar-On]

Hi Alan!

I'm using the code that you helped me with regarding the PB equation in this equation of rho. I solved it correctly for one set of kdim and Bdim but for some reason that I don't understand It doesn't run as a loop for a array of these parameters. This is my code:

tolerance = 10^-3;  % Set as desired
kdimarray = 0.5:0.5:3; Bdimarray = 0.5:0.5:3;
xspan = [0 3];  % inf = 5 here!
% Initial guess for gradients  dpsidx(x=0)
%dp1 = -1;
dp1 = [-1,-1,-1,-1,-1,-1];
delta = -10^-10;  
flag = 1; % Repeat indicator
itmax = 100; % Maximum allowable number of iterations
its = 0;
for i = 1:6
    kdim = kdimarray(i); Bdim = Bdimarray(i);
while (flag==1) && its<itmax 
    % Call ode solver twice with different guesses for initial gradient
    % and get values of psi at boundary in return psi(x=inf)
    [~, p1] = BP([0.74 dp1(i)],xspan);
    P1 = p1(end);
    [x, p2] = BP([0.74 (dp1(i)+delta)],xspan);
    P2 = p2(end);
       
    if abs(P2)<tolerance  % then we are finished
          flag = 0;
    else  % we need another iteration
          % Secant improvement
          dP = P2 - P1;  % Difference between end values of psi
          dp1 = dp1 - delta*P1/dP; 


    end
    
    its = its + 1;
end
end
% Display data
format long
disp('rho(x=inf) and drhodx(x=0)')
disp([P2 dp1])
disp(['number of iterations = ' num2str(its)])


 plot(x,p2(:,1),'b')
 hold on
 xlabel('x')
 ylabel('\rho')


function [x,p] = BP(y0,xspan)
%Shooting method for the BP eq
options = odeset('RelTol',1e-3,'AbsTol',1e-3);
kdimarray = 0.5:0.5:3; Bdimarray = 0.5:0.5:3;
 kdim = kdimarray(i); Bdim = Bdimarray(i);
[x,p] = ode45(@eqs, xspan, y0 , kdim, Bdim, [], options);

%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2)),grid
%legend('\psi','d\psi dx')
end

function dy = eqs(~,y,kdim,Bdim)  % x not used within the function so replaced by ~
%PB equation
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim./kdim;
end

Can you help? I'm pretty sure it is something tiny that I've missed.

Thanks,

Roi

[Alan Stevens]

Well, the following works, but only you can tell if the resuts are sensible.

tolerance = 10^-3;  % Set as desired
kdimarray = 0.5:0.5:3; Bdimarray = 0.5:0.5:3;
xspan = linspace(0,3,50);  %%%%%  To get same number of points for each iteration
% Initial guess for gradients  dpsidx(x=0)
%dp1 = -1;
dp1 = [-1,-1,-1,-1,-1,-1];
delta = -10^-10;  

itmax = 100; % Maximum allowable number of iterations

for i = 1:numel(kdimarray)
    flag = 1;             %%%%%%%%%% Need this inside the for loop
    its = 0;
    kdim = kdimarray(i); Bdim = Bdimarray(i);
    while (flag==1) && its<itmax 
        % Call ode solver twice with different guesses for initial gradient
        % and get values of psi at boundary in return psi(x=inf)
        [~, p1] = BP([0.74 dp1(i)],xspan,kdim,Bdim);
        P1 = p1(end);
        [x, p2] = BP([0.74 (dp1(i)+delta)],xspan,kdim,Bdim);
        P2 = p2(end);

        if abs(P2)<tolerance  % then we are finished
              flag = 0;
        else  % we need another iteration
              % Secant improvement
              dP = P2 - P1;  % Difference between end values of psi
              dp1(i) = dp1(i) - delta*P1/dP; 
        end

        its = its + 1;
    end
Psol(i) = P2;         %%%%%%%%%%  Save values seperately for each loop 
Dp1(i) = dp1(i);      %%%%%%%%%%       Ditto
rho(:,i) = p2(:,1);   %%%%%%%%%%       Ditto

end
% Display data
format long
disp('rho(x=inf) and drhodx(x=0)')
disp([Psol;Dp1])

 plot(x,rho),grid
 hold on
 xlabel('x')
 ylabel('\rho')
 legend('i = 1','i = 2','i = 3','i = 4','i = 5','i = 6')


function [x,p] = BP(y0,xspan,kdim,Bdim)
%Shooting method for the BP eq
options = odeset('RelTol',1e-3,'AbsTol',1e-3);
[x,p] = ode45(@(x,y)eqs(xspan, y0,kdim,Bdim), xspan, y0 ,options);  % Pass extra parameters to ode45

%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2)),grid
%legend('\psi','d\psi dx')
end

function dy = eqs(~,y,kdim,Bdim)  % x not used within the function so replaced by ~
%PB equation  
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim./kdim;
end

[Roi Bar-On]

Hey Alan, Thanks again for the help.

There's something weird.. If I run my code without any loop for a single set of kdim and Bdim I'm getting physical results and a smooth curve. In this case I'm getting linear curves with some negative values for rho (not physical).

This is the code that I'm using for a singke computation:

tolerance = 10^-3;  % Set as desired
xspan = [0 3];  % inf = 5 here!
% Initial guess for gradients  dpsidx(x=0)
dp1 = -1;
delta = -10^-10;  
flag = 1; % Repeat indicator
itmax = 100; % Maximum allowable number of iterations
its = 0;
while (flag==1) && its<itmax
    % Call ode solver twice with different guesses for initial gradient
    % and get values of psi at boundary in return psi(x=inf)
    [~, p1] = BP([0.74 dp1],xspan);
    P1 = p1(end);
    [x, p2] = BP([0.74 (dp1+delta)],xspan);
    P2 = p2(end);
       
    if abs(P2)<tolerance  % then we are finished
          flag = 0;
    else  % we need another iteration
          % Secant improvement
          dP = P2 - P1;  % Difference between end values of psi
          dp1 = dp1 - delta*P1/dP;


    end
   
    its = its + 1;
end
% Display data
format long
disp('rho(x=inf) and drhodx(x=0)')
disp([P2 dp1])
disp(['number of iterations = ' num2str(its)])


 plot(x,p2(:,1),'b')
 xlabel('x'),ylabel('\rho')


function [x,p] = BP(y0,xspan)
%Shooting method for the BP eq
options = odeset('RelTol',1e-3,'AbsTol',1e-3);
[x,p] = ode45(@eqs, xspan, y0 ,options);

%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2)),grid
%legend('\psi','d\psi dx')
end

function dy = eqs(~,y)  % x not used within the function so replaced by ~
%PB equation
kdim=2;Bdim=2;
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim./kdim;
end

I don't understand what happend. for kdim=2 and Bdim=2 the curves are different between the two codes..

[Alan Stevens]

I think your log term in function eqs is causing problems. Although your code above appears to work (giving a curved plot) the values it produces for rho are complex! (The plot ignores the imaginary components). This can be overcome by setting dp1 = -0.1 instead of -1. However, you still max out on iterations!

Will look further to see why the other code is weird.

[Roi Bar-On]

In the ''regular'' code it's only 6 interations

[Alan Stevens]

What is the mathematical form of your current 2nd order ode?

[Roi Bar-On]

By the way, how did you know that -0.1 will not produce complex numbers? Which kind of calculations have you done?

[Roi Bar-On]

by the way 0.25 can be any number below 0.74 (we still haven't determined it at this moment..)

[Alan Stevens]

Are you sure? Your original had B as the last term not B/K, and the first term wasn't multiplied by rho.

[Roi Bar-On]

Yes, in the original post I had a mistake which I have fixed in previous comments.

[Alan Stevens]

OK. So here's something to think about

So, which is more likely, B positive or rho(inf) less than 0.5?

[Roi Bar-On]

Wow.. Amazing job. Basically we don't the value of rho_inf but it can be below 0.5.. Also, B can be negative (in theory) but in our work for now it is positive. So rho_inf can be found between 0.5 to 0.74.

Alain, how could you guess that initial gradient of -0.1 will avoid complex numbers?

[Alan Stevens]

"how could you guess that initial gradient of -0.1 will avoid complex numbers?"

Nothing clever - I just tried it!

[Alan Stevens]

I notice that if B and K have opposite signs it looks like there is an oscillating solution, which would mean there is no constant limit at infinity for rho.

[Roi Bar-On]

Have you seen it graphically?

[Alan Stevens]

This

B = 2; K = -2;
rho0 = 0.74;
drhodx0 = 0;
xspan = 0:0.01:5;
IC = [rho0 drhodx0];

[x, R] = ode45(@(x,R) odefn(x,R,B,K),xspan,IC);
rho = R(:,1);
drhodx = R(:,2);

plot(x,rho,x,drhodx),grid
xlabel('x'),ylabel('\rho & d\rho dx')
legend('\rho','d\rho dx')

function dRdx = odefn(~,R,B,K)       
         rho    = R(1);
         drhodx = R(2);
         dRdx = [drhodx;
                (2*B*rho + log(rho) - B)/K];
end

produces this

[Roi Bar-On]

Okay I see. So if I work now with large enough values of B2 and K , the code (with a for loop with arrays of B2 and K) should yield decent plots? Or I have to choose carefully my initial gradient each time?

[Roi Bar-On]

Because it doesn't work well..

[Alan Stevens]

My previous post wasn't about large values of B and K, but about the fact that you can get oscillatory behaviour if they are of opposite sign.

I can't advise on what values to use as I don't know the background to the problem (Is it the model of a physical or chemical system, for example, or just something of mathematical interest?).

You still need to choose a value for the initial gradient.

[Roi Bar-On]

That's not what I'm talking about.. Regardless of the values i;ve asked:

There's something weird.. If I run my code without any loop for a single set of kdim and Bdim I'm getting physical results and a smooth curve. In this case I'm getting linear curves with some negative values for rho (not physical).

This is the code that I'm using for a single computation:

tolerance = 10^-3;  % Set as desired
xspan = [0 3];  % inf = 5 here!
% Initial guess for gradients  dpsidx(x=0)
dp1 = -1;
delta = -10^-10;  
flag = 1; % Repeat indicator
itmax = 100; % Maximum allowable number of iterations
its = 0;
while (flag==1) && its<itmax
    % Call ode solver twice with different guesses for initial gradient
    % and get values of psi at boundary in return psi(x=inf)
    [~, p1] = BP([0.74 dp1],xspan);
    P1 = p1(end);
    [x, p2] = BP([0.74 (dp1+delta)],xspan);
    P2 = p2(end);
       
    if abs(P2)<tolerance  % then we are finished
          flag = 0;
    else  % we need another iteration
          % Secant improvement
          dP = P2 - P1;  % Difference between end values of psi
          dp1 = dp1 - delta*P1/dP;
    end
   
    its = its + 1;
end
% Display data
format long
disp('rho(x=inf) and drhodx(x=0)')
disp([P2 dp1])
disp(['number of iterations = ' num2str(its)])
 plot(x,p2(:,1),'b')
 xlabel('x'),ylabel('\rho')
function [x,p] = BP(y0,xspan)
%Shooting method for the BP eq
options = odeset('RelTol',1e-3,'AbsTol',1e-3);
[x,p] = ode45(@eqs, xspan, y0 ,options);
%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2)),grid
%legend('\psi','d\psi dx')
end
function dy = eqs(~,y)  % x not used within the function so replaced by ~
%PB equation
kdim=2;Bdim=2;
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim./kdim;
end

I don't understand what happend. for kdim=2 and Bdim=2 the curves are different between the two codes..

I just want the code to work well with the for loop and profuce same plots as I get in the ''regular'' code

[Alan Stevens]

Try this

tolerance = 10^-3;  % Set as desired
xspan = 0:0.1:3;  % inf = 5 here!
kdim = [2 2.1 2.2]; Bdim = [2 1.9 1.8];                    %**************%
rho = zeros(numel(xspan),numel(kdim));                     %**************%
drhodx = zeros(numel(xspan),numel(kdim));                  %**************%
for i = 1:numel(kdim)                                      %**************%
    kd = kdim(i); Bd = Bdim(i);                            %**************%
    % Initial guess for gradients  dpsidx(x=0)
    dp1 = -0.1;   % This value won't work for all combinations of k & B
    delta = -10^-10;  
    flag = 1; % Repeat indicator
    itmax = 100; % Maximum allowable number of iterations
    its = 0;
    while (flag==1) && its<itmax
        % Call ode solver twice with different guesses for initial gradient
        % and get values of psi at boundary in return psi(x=inf)
        [~, p1] = BP([0.74 dp1],xspan,kd,Bd);              %**************%
        P1 = p1(end,2);
        [x, p2] = BP([0.74 (dp1+delta)],xspan,kd,Bd);      %**************%
        P2 = p2(end,2);

        if abs(P2)<tolerance  % then we are finished
              flag = 0;
        else  % we need another iteration
              % Secant improvement
              dP = P2 - P1;  % Difference between end values of psi
              dp1 = dp1 - delta*P1/dP;
        end

        its = its + 1;
    end
    rho(:,i) = p2(:,1);
    drhodx(:,i) = p2(:,1);
end
% Display data
% format long
% disp('rho(x=inf) and drhodx(x=0)')
% disp([P2 dp1])
% disp(['number of iterations = ' num2str(its)])
 plot(x,rho),grid
 xlabel('x'),ylabel('\rho')
 
function [x,p] = BP(y0,xspan,kd,Bd)
%Shooting method for the BP eq
options = odeset('RelTol',1e-3,'AbsTol',1e-3);
[x,p] = ode45(@(x,y)eqs(x,y,kd,Bd), xspan, y0 ,options);      %***********%
%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2)),grid
%legend('\psi','d\psi dx')
end
function dy = eqs(~,y,kd,Bd)                               %**************%
%PB equation
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 2.*Bd./kd.*y(1)+kd.^(-1).*log(y(1))-Bd./kd;
end

However, note that although your single value code seemed to work ok, it generated complex values that were ignored by the plot command. The above works ok, with no complex numbers. However, not all combinations will do so with the single value used for the single initial value of dpi. You will probably have to change it for some combinations of k and B.

[Alan Stevens]

Replacing the eqs function with the following

function dy = eqs(~,y,kd,Bd)                               
%PB equation
dy = zeros(2,1);    % a column vector
lny = log(abs(y(1)));
dy(1) = y(2);
dy(2) = 2.*Bd./kd.*y(1)+kd.^(-1).*lny-Bd./kd;
end

removes the complex number problem.

[Alan Stevens]

Also, it's a good idea to insert

    p1 = zeros(numel(xspan),2);
    p2 = zeros(numel(xspan),2);

immediately before the

while (flag==1) && its<itmax

line.

[Roi Bar-On]

Hi Alain.

You guys helped me a lot with my previous code. I've implemented some changes and improvements to it but now having some new problems. This is still my ODE with the following bc's:

My code includes the shooting method and an intellectual guess of the derivative of rho at x=0 (transforming one BVP problems to one IVP problem), using the Implicit Euler method to write down the derivatives, solve the algebraic system with the Newton-Raphson method and make sure that I converge with infinity norm condition. I try to optimize my results by conducting bisection on the difference between a current value of rho to the its value at the bc that I chose (minimizing f1). This is the implicit Euler code:

function[x,y,f1,BISEC]=implicit_euler_method_DFT(y10,y20,x0,xf,dx,m,n)%i.g=-0.4394
x=x0:dx:xf;
y=zeros(2,length(x));
y(1,1)=y10;
y(2,1)=y20;
B=1;k=1;bc=0.25;
for i=1:1:length(x)-1
    y(:,i+1)=y(:,i);
    J=Jacovian(y(:,i));
    g=[y(1,i+1)-dx*y(2,i+1)-y(1,i);
        y(2,i+1)-dx*(2*B/k*y(1,i+1)+1/k*log(abs(y(1,i+1)))-B/k)-y(2,i)];
    d=gausselemination(J,g);
    f1=y(1,i+1)-bc;
    BISEC=bisectionDFTbc(f1);
    while (infinity_norm(d)>10^(-m))||(infinity_norm(g)>10^(-n))
        y(:,i+1)=y(:,i+1)+d;
        J=Jacovian(y(:,i+1));
        g=[y(1,i+1)-dx*y(2,i+1)-y(1,i);
        y(2,i+1)-dx*(2*B/k*y(1,i+1)+1/k*log(abs(y(1,i+1)))-B/k)-y(2,i)];
        d=gausselemination(J,g);
        f1=y(1,i+1)-bc;
        BISEC=bisectionDFTbc(f1);
    end
end
end

where x is the x tilde basically, y is rho tilde, y10 is the initial value for rho, y20 is the initial guess (-0.4394) for its derivative at x=0,x0=0, xf can be 1 or more, dx is to our desire, m and n are tolerances. The 2nd code for the bisection is:

function [x,iter_num]= bisectionDFTbc(a,b,n,m)
 iter_num=0;%iteration number
 x=(a+b)/2;
 while sqrt((a-b)^2)>10^(-m)||sqrt(f1(x)^2)>10^(-n)
     iter_num=iter_num+1;
     x=(a+b)/2;
     if f1(a)*f1(x)<0
         b=x;
     else
         a=x;
     end
 end
 x=(b+a)/2;
end

I have two main problems:

1.The y values go to infinity very fast - I don't get a stable solution.

2.Obviously I don't connect the two functions well.

I would appreciate your help.

Thanks,

Roi

Answer 2

[Roi Bar-On]

That's what I've done before. I've putted it back.

Thank you!