Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss
Questions regarding bvp4c with unknown parameters
188 views
Skip to first unread message
Lauren
Mar 14, 2012, 9:50:17 PM3/14/12
to
Hello,
I am trying to learn how to use the bvp4c solver to solve a second order differential equation. I've been around the documentation and prior newsgroup postings and I really just keep getting confused since the examples seem to be either to complex or don't make sense to me.
I'm using a simplified version of the equation I ultimately want to solve so that I can actually learn how to do this.
The 2nd order ODE is C" = p(1) * exp(-p(2) * x) - p(3)*exp(-p(4)*x)*C where p(1:4) are unknown parameters
The boundary conditions are that C(0) = Cbw and C(xm) = G where xm is the last value in a matrix of data and G is a scalar computed from data.
I've rewritten the 2nd order ODE as:
y(1) = C
y(1)' = C' = y(2)
y(2) = C" = p(1) * exp(-p(2) * x) - p(3)*exp(-p(4)*x)*y(1)
and my boundary conditions are written as:
ya(1) - Cbw
ya(2) - Gm
My questions regarding this process are:
1) Did I actually write the system of equations and boundary conditions correctly? (I really don't understand how the examples and tutorials write boundary conditions)
2) I am getting an error that says " Attempted to access y(2); index out of bounds because numel(y)=1" and I am not sure why. I've tried to research this error, but I haven't been able to figure out how to apply the suggestions from prior posts to my code
3) Once the ODE is solved, how do I go about optimizing the parameters to the set of data?
Thank you.
~~~~~~~~~~~~~~~~~
function numeric
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
%define the concentration gradient at deepest depth
global Gm
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
%initial guess for parameters
p = [.000005 .000000005];
solinit = bvpinit(linspace(0,50,250),0.022,p);
sol = bvp4c(@ode, @bc, solinit);
function dydx = ode(x,y,p)
dydx= [y(2);p(1)*exp(-p(2)*x) - p(3)*exp(-p(4)*x)*y(1)];
function boundary = bc(ya,yb,p)
N0=0.022;
global Gm
boundary = [ya(1) - N0
yb(2)- Gm];
I am trying to learn how to use the bvp4c solver to solve a second order differential equation. I've been around the documentation and prior newsgroup postings and I really just keep getting confused since the examples seem to be either to complex or don't make sense to me.
I'm using a simplified version of the equation I ultimately want to solve so that I can actually learn how to do this.
The 2nd order ODE is C" = p(1) * exp(-p(2) * x) - p(3)*exp(-p(4)*x)*C where p(1:4) are unknown parameters
The boundary conditions are that C(0) = Cbw and C(xm) = G where xm is the last value in a matrix of data and G is a scalar computed from data.
I've rewritten the 2nd order ODE as:
y(1) = C
y(1)' = C' = y(2)
y(2) = C" = p(1) * exp(-p(2) * x) - p(3)*exp(-p(4)*x)*y(1)
and my boundary conditions are written as:
ya(1) - Cbw
ya(2) - Gm
My questions regarding this process are:
1) Did I actually write the system of equations and boundary conditions correctly? (I really don't understand how the examples and tutorials write boundary conditions)
2) I am getting an error that says " Attempted to access y(2); index out of bounds because numel(y)=1" and I am not sure why. I've tried to research this error, but I haven't been able to figure out how to apply the suggestions from prior posts to my code
3) Once the ODE is solved, how do I go about optimizing the parameters to the set of data?
Thank you.
~~~~~~~~~~~~~~~~~
function numeric
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
%define the concentration gradient at deepest depth
global Gm
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
%initial guess for parameters
p = [.000005 .000000005];
solinit = bvpinit(linspace(0,50,250),0.022,p);
sol = bvp4c(@ode, @bc, solinit);
function dydx = ode(x,y,p)
dydx= [y(2);p(1)*exp(-p(2)*x) - p(3)*exp(-p(4)*x)*y(1)];
function boundary = bc(ya,yb,p)
N0=0.022;
global Gm
boundary = [ya(1) - N0
yb(2)- Gm];
Torsten
Mar 15, 2012, 6:03:01 AM3/15/12
to Lauren
> solinit = bvpinit(linspace(0,50,250),0.022,p);
You supply only one value (0.022) for the solution y.
But y has 2 components:
solinit = bvpinit(linspace(0,50,250),[0.022 ??]);
For ?? you will have to supply a guess for C'.
> sol = bvp4c(@ode, @bc, solinit);
>
sol = bvp4c(@(x,y)ode(x,y,p), @(ya,yb)bc(ya,yb,p), solinit);
> function dydx = ode(x,y,p)
> dydx= [y(2);p(1)*exp(-p(2)*x) - p(3)*exp(-p(4)*x)*y(1)];
>
> function boundary = bc(ya,yb,p)
> N0=0.022;
> global Gm
> boundary = [ya(1) - N0
> yb(2)- Gm];
the last condition must read
boundary = [ya(1) - N0
yb(1)- Gm];
First try to run bvp4c with a fixed p-vector.
If this works fine, you will have to couple
bvp4c and a least-squares solver (e.g. lsqcurvefit)
for parameter fitting.
Best wishes
Torsten.
Lauren
Mar 15, 2012, 12:28:13 PM3/15/12
to
Hi Torsten,
Thank you so much for your response.
For your first question about the number of parameters - that was a copy-paste error, I had been trying to see if simplifying the equation further would help (it didn't).
I have implemented your suggestions, but now it says that there are too many input arguments in this line:
sol = bvp4c(@(x,y)ode, @(ya,yb)bc(ya,yb,p), solinit);
Where the error is:
??? Error using ==> numeric>@(x,y)ode
Too many input arguments.
When I change it back to @ode it still says there are too many input arguments:
??? Error using ==> numeric>@(ya,yb)bc(ya,yb,p)
Too many input arguments.
Also, in the line: solinit = bvpinit(linspace(0,50,250),[0.022,0.05],p);
In your suggestion I noticed you left out the 'p' term, but if I do that, it says that y is undefined, so I put it back in.
Finally, with regards to the boundary condition function, why is the second part yb(1) -Gm instead of yb(2)-Gm? I understand that yb is indicating y(b) or y at the rightmost boundary, but how do you choose (1) or (2)?
My updated code is pasted below.
Thank you!
p = [.000005 .000000005 .000005 .000000005]; %initial guess for parameters p1-4
solinit = bvpinit(linspace(0,50,250),[0.022,0.05],p);
sol = bvp4c(@(x,y)ode, @(ya,yb)bc(ya,yb,p), solinit);
function dydx = ode(x,y,p)
dydx= [y(2);p(1)*exp(-p(2)*x) - p(3)*exp(-p(4)*x)*y(1)];
function boundary = bc(ya,yb,p)
N0=0.022;
global Gm
boundary = [ya(1) - N0
yb(1)- Gm];
Thank you so much for your response.
For your first question about the number of parameters - that was a copy-paste error, I had been trying to see if simplifying the equation further would help (it didn't).
I have implemented your suggestions, but now it says that there are too many input arguments in this line:
sol = bvp4c(@(x,y)ode, @(ya,yb)bc(ya,yb,p), solinit);
Where the error is:
??? Error using ==> numeric>@(x,y)ode
Too many input arguments.
When I change it back to @ode it still says there are too many input arguments:
??? Error using ==> numeric>@(ya,yb)bc(ya,yb,p)
Too many input arguments.
Also, in the line: solinit = bvpinit(linspace(0,50,250),[0.022,0.05],p);
In your suggestion I noticed you left out the 'p' term, but if I do that, it says that y is undefined, so I put it back in.
Finally, with regards to the boundary condition function, why is the second part yb(1) -Gm instead of yb(2)-Gm? I understand that yb is indicating y(b) or y at the rightmost boundary, but how do you choose (1) or (2)?
My updated code is pasted below.
Thank you!
solinit = bvpinit(linspace(0,50,250),[0.022,0.05],p);
sol = bvp4c(@(x,y)ode, @(ya,yb)bc(ya,yb,p), solinit);
function dydx = ode(x,y,p)
dydx= [y(2);p(1)*exp(-p(2)*x) - p(3)*exp(-p(4)*x)*y(1)];
function boundary = bc(ya,yb,p)
N0=0.022;
global Gm
boundary = [ya(1) - N0
Torsten
Mar 16, 2012, 3:54:53 AM3/16/12
to Lauren
try this:
function numeric
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
%define the concentration gradient at deepest depth
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
N0=0.022;
p = [.000005 .000000005 .000005 .000000005]; %initial guess for parameters p1-4
solinit = bvpinit(linspace(0,50,250),[0.022,0.05]);
sol = bvp4c(@ode, @bc, solinit);
function dydx = ode(x,y)
global Gm N0 p
dydx= [y(2);p(1)*exp(-p(2)*x) - p(3)*exp(-p(4)*x)*y(1)];
function boundary = bc(ya,yb)
global Gm N0 p
boundary = [ya(1) - N0
yb(1)- Gm];
Concerning your last question:
yb(1)- Gm];
There is contradictory information in your first post.
If you want to set C(xm) = G (as you write in your introduction),
you will have to set yb(1)-Gm.
If you want to set dC/dx = Gm (as I conclude from the definition of Gm),
you will have to set yb(2)-Gm.
Best wishes
Torsten.
Lauren
Mar 16, 2012, 6:57:12 PM3/16/12
to
Hi Torsten,
Thank you very much! The new changes to the code appear to work, and I do finally understand where the boundary conditions come from. (You are correct that Gm is equivalent to dC/dx at the right hand boundary).
Now, regarding the 3rd question from my original post, I am attempting to optimize the parameters in 'p' to fit the dataset. I can think of two ways of doing this, and neither one seems to work.
First method:
My understanding of bvp4c from the documentation is that it can find 'unknown parameters p of the form y' = f(x,y,p)' Using the command 'sol.p' to have it tell me what the 'final' parameters are doesn't appear to work, so it appears that bvp4c is treating p as a vector of constants rather than as an initial guess for solving the ODE since p does not change after running the solver.
When I try adding p to the line for solinit to define p as a vector of parameters, it gives me an error saying there are too many input parameters.
One example I found seems to indicate that if there are is an unknown parameter, more boundary conditions are required. I have some more of these, but they apply to a boundary between a and b. (xhat)
Second method:
In theory, I think I can just use p as a set of 'constants' for solving the equation, and then optimize them in lsqcurvefit, but I'm not really sure how to do that. I've used lsqcurvefit when starting from the analytical solution to an ODE, but I am not sure that I'm passing lsqcurvefit the correct function.
The code pasted below includes the full function I am trying to solve (which bvp4c is able to solve) along with my attempts at getting lsqcurvefit to work.
The error it is giving me is that "function value and ydata sizes are incommensurate" and thinking about it, I can see exactly why that isn't working (y and y data have different sizes), but I am not sure how to change the input to lsqcurvefit to give it the correct function to fit to.
Code is below, Thanks.
~~~~~~~~~~~~~~~~~~~~~~
ydata = Data(:,2);
%define the constants, etc
global Gm N0 xhat p
%concentration gradient at deepest depth
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
%concentration at the surface
N0=0.022;
%boundary between production and consumption
xhat = 18.4;
%vector of initial guesses for the parameters
p = [.0000076 .00000004 .000002 .000000006 .000005 .000000005]; %initial guess for parameters p1-6
solinit = bvpinit(linspace(0,50,250),[0.022,0.03]);
figure(1)
plot(sol.x,sol.y,'b-');
hold on
plot(xdata,ydata,'go');
xlabel('Depth')
ylabel('Nitrate (µol/cm3)')
legend('Nitrate model','Nitrate data')
%Optimize Parameters in P to fit the data
options=optimset('FunValCheck','on','MaxIter',1000,'MaxFunEvals',100000,'TolX',1e-14,'Tolfun',1e-14);
pfinal=lsqcurvefit(@ode,p,xdata,ydata,[],[],options);
%The ODE
function dydx = ode(x,y)
global Gm N0 xhat p
f=zeros(length(x),1);
for i=1:length(x)
if (x(i) <= xhat)
A = 1;
else
A = 0;
end
end
dydx= [y(2); A*(p(1)*exp(-p(2)*x(i))+ p(3)*exp(-p(4)*x(i)) ) - (1-A)*(p(5)*exp(-p(6)*x(i)))*y(1)];
%The boundary conditions
function boundary = bc(ya,yb)
global Gm N0 xhat p
Thank you very much! The new changes to the code appear to work, and I do finally understand where the boundary conditions come from. (You are correct that Gm is equivalent to dC/dx at the right hand boundary).
Now, regarding the 3rd question from my original post, I am attempting to optimize the parameters in 'p' to fit the dataset. I can think of two ways of doing this, and neither one seems to work.
First method:
My understanding of bvp4c from the documentation is that it can find 'unknown parameters p of the form y' = f(x,y,p)' Using the command 'sol.p' to have it tell me what the 'final' parameters are doesn't appear to work, so it appears that bvp4c is treating p as a vector of constants rather than as an initial guess for solving the ODE since p does not change after running the solver.
When I try adding p to the line for solinit to define p as a vector of parameters, it gives me an error saying there are too many input parameters.
One example I found seems to indicate that if there are is an unknown parameter, more boundary conditions are required. I have some more of these, but they apply to a boundary between a and b. (xhat)
Second method:
In theory, I think I can just use p as a set of 'constants' for solving the equation, and then optimize them in lsqcurvefit, but I'm not really sure how to do that. I've used lsqcurvefit when starting from the analytical solution to an ODE, but I am not sure that I'm passing lsqcurvefit the correct function.
The code pasted below includes the full function I am trying to solve (which bvp4c is able to solve) along with my attempts at getting lsqcurvefit to work.
The error it is giving me is that "function value and ydata sizes are incommensurate" and thinking about it, I can see exactly why that isn't working (y and y data have different sizes), but I am not sure how to change the input to lsqcurvefit to give it the correct function to fit to.
Code is below, Thanks.
~~~~~~~~~~~~~~~~~~~~~~
function numeric
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
xdata = Data(:,1);
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
ydata = Data(:,2);
%define the constants, etc
global Gm N0 xhat p
%concentration gradient at deepest depth
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
%concentration at the surface
N0=0.022;
%boundary between production and consumption
xhat = 18.4;
%vector of initial guesses for the parameters
p = [.0000076 .00000004 .000002 .000000006 .000005 .000000005]; %initial guess for parameters p1-6
solinit = bvpinit(linspace(0,50,250),[0.022,0.03]);
sol = bvp4c(@ode, @bc, solinit);
%plot the solution to the ODE from bvp4c
figure(1)
plot(sol.x,sol.y,'b-');
hold on
plot(xdata,ydata,'go');
xlabel('Depth')
ylabel('Nitrate (µol/cm3)')
legend('Nitrate model','Nitrate data')
%Optimize Parameters in P to fit the data
options=optimset('FunValCheck','on','MaxIter',1000,'MaxFunEvals',100000,'TolX',1e-14,'Tolfun',1e-14);
pfinal=lsqcurvefit(@ode,p,xdata,ydata,[],[],options);
%The ODE
function dydx = ode(x,y)
global Gm N0 xhat p
f=zeros(length(x),1);
for i=1:length(x)
if (x(i) <= xhat)
A = 1;
else
A = 0;
end
end
dydx= [y(2); A*(p(1)*exp(-p(2)*x(i))+ p(3)*exp(-p(4)*x(i)) ) - (1-A)*(p(5)*exp(-p(6)*x(i)))*y(1)];
%The boundary conditions
function boundary = bc(ya,yb)
global Gm N0 xhat p
boundary = [ya(1) - N0
yb(2)- Gm];
Torsten
Mar 17, 2012, 7:53:23 AM3/17/12
to Lauren
I'll show you below the rough way - you will have to fill in the details.
options=optimset('FunValCheck','on','MaxIter',1000,'MaxFunEvals',100000,'TolX',1e-14,'Tolfun',1e-14);
pfinal=lsqcurvefit(@myfun,p,xdata,ydata,[],[],options)
function F = myfun(x,xdata)
global Gm N0 xhat p
p=x;
solinit = bvpinit(linspace(0,50,250),[0.022,0.03]);
sol = bvp4c(@ode, @bc, solinit);
% There is a function in bvp4c that interpolates the solution
sol = bvp4c(@ode, @bc, solinit);
% sol to your data points xdata. Use this function here to get
% the array ydata_simulated.
ydata_simulated = ...;
F = ydata_simulated;
%The ODE
function dydx = ode(x,y)
global Gm N0 xhat p
f=zeros(length(x),1);
for i=1:length(x)
if (x(i) <= xhat)
A = 1;
else
A = 0;
end
end
%The boundary conditions
function boundary = bc(ya,yb)
global Gm N0 xhat p
boundary = [ya(1) - N0
yb(2)- Gm];
Best wishes
Torsten.
Lauren
Mar 21, 2012, 2:54:33 PM3/21/12
to
Hi Torsten,
Just wanted to report that I was able to figure out how to get the lsqcurvefit-bvp4c connection to work!
Thank you so much for your help!
~ Lauren
Just wanted to report that I was able to figure out how to get the lsqcurvefit-bvp4c connection to work!
Thank you so much for your help!
~ Lauren
Lauren
Mar 22, 2012, 6:34:31 PM3/22/12
to
Hello again,
So everything seemed to be working great, and then I realized that I left out an important constant from the ODE. When I added this (1/Dno3) to the code, now it keeps running and never seems to converge on a 'final' solution.
I've tried turning on the options in bvpset so that I can see the stats as it's running and change the tolerance, but it never comes to a solution (I let it run overnight)
I've tried reducing the number of points it generates, changing the initial guess, and reducing the tolerance, but nothing seems to get it to converge anymore and it appears to be stuck in some sort of loop. The bvp4c stats keep coming back with nearly the same numbers each time for mesh points, residual, and number of calls.
Any advice on how to get the solver 'unstuck' would be much appreciated.
My code is below, but the only thing that has really changed from the last posting is adding (1/Dno3) to the ODE function
dydx= [y(2); (1/Dno3)*((A)*(p(1,1)*exp(-p(1,2)*x(i))+ p(1,3)*exp(-p(1,4)*x(i)) ) - (1-A)*(p(1,5)*exp(-p(1,6)*x(i)))*y(1))];
~~~~~
function numeric
%define the constants, etc
global Gm N0 xhat p xmin xmax g Dno3
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
xdata = Data(:,1);
ydata = Data(:,2);
xmin=min(Data(:,1));
xmax=max(Data(:,1));
%concentration gradient at deepest depth
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
Dno3=4.2e-06;
DO2=5.1e-06;
O2zero=.240;
N0=0.022; %NO3 concentration at the surface
g = .1; %gamma, the ratio of NO3/O2
xhat = 18.4;%boundary between production and consumption
%vector of initial guesses for the parameters
p = [.00076 .00004 .0002 .00006 .0005 .00005]; %initial guess for parameters p1-6
%Setup to fit the data with the solution of the ODE
options=optimset('FunValCheck','on','MaxIter',5000,'MaxFunEvals',100000,'TolX',1e-8,'Tolfun',1e-8);
[pfinal,resnorm,exitflag,output]=lsqcurvefit(@fitN,p,xdata,ydata,[],[],options);
%Plot output of lsqcurvefit
xsim = linspace(xmin,xmax,50);
yout = fitN(pfinal,xsim);
exitflag
output
pfinal
%plot the solution to the ODE from bvp4c
figure(1)
plot(xsim,yout,'b-');
function F = fitN(x,xdata)
global Gm N0 xhat p xmin xmax Dno3
p=x;
options= bvpset('Stats','on','RelTol',1e-4,'abstol',1e-3);
solinit = bvpinit(linspace(xmin,xmax,50),[0.022,0.03]);
sol = bvp4c(@ode, @bc, solinit,options);
ydata_sim = deval(sol,xdata)'; %evaluates solution of bvp at the xdata points
F = ydata_sim(:,1);
%The ODE for Nitrogen
function dydx = ode(x,y)
global Gm N0 xhat p xmin xmax Dno3
f=zeros(length(x),1);
for i=1:length(x)
if (x(i) <= xhat)
A = 1;
else
A = 0;
end
end
dydx= [y(2); (1/Dno3)*((A)*(p(1,1)*exp(-p(1,2)*x(i))+ p(1,3)*exp(-p(1,4)*x(i)) ) - (1-A)*(p(1,5)*exp(-p(1,6)*x(i)))*y(1))];
%The boundary conditions for Nitrogen
function boundary = bc(ya,yb)
global Gm N0 xhat p xmin xmax Dno3
So everything seemed to be working great, and then I realized that I left out an important constant from the ODE. When I added this (1/Dno3) to the code, now it keeps running and never seems to converge on a 'final' solution.
I've tried turning on the options in bvpset so that I can see the stats as it's running and change the tolerance, but it never comes to a solution (I let it run overnight)
I've tried reducing the number of points it generates, changing the initial guess, and reducing the tolerance, but nothing seems to get it to converge anymore and it appears to be stuck in some sort of loop. The bvp4c stats keep coming back with nearly the same numbers each time for mesh points, residual, and number of calls.
Any advice on how to get the solver 'unstuck' would be much appreciated.
My code is below, but the only thing that has really changed from the last posting is adding (1/Dno3) to the ODE function
dydx= [y(2); (1/Dno3)*((A)*(p(1,1)*exp(-p(1,2)*x(i))+ p(1,3)*exp(-p(1,4)*x(i)) ) - (1-A)*(p(1,5)*exp(-p(1,6)*x(i)))*y(1))];
~~~~~
function numeric
%define the constants, etc
global Gm N0 xhat p xmin xmax g Dno3
%Data - column 1 is depth, column 2 is concentration in umol/cm3
Data = ...
[0 0.0220
0.12 0.0235
0.28 0.0248
0.47 0.0275
0.65 0.0287
0.82 0.0289
1.02 0.0292
1.21 0.0295
1.905 0.0278
3.81 0.0310
5.72 0.0335
11.42 0.0310
13.32 0.0267
15.23 0.0225
18.08 0.0180
22.85 0.0096
27.61 0.0043
27.61 0.0043
32.37 0.0007
41.90 0.0004];
xdata = Data(:,1);
ydata = Data(:,2);
xmax=max(Data(:,1));
%concentration gradient at deepest depth
Gm = (Data(end-1,2)-Data(end,2))/(Data(end-1,1)-Data(end,1));
DO2=5.1e-06;
O2zero=.240;
N0=0.022; %NO3 concentration at the surface
g = .1; %gamma, the ratio of NO3/O2
xhat = 18.4;%boundary between production and consumption
%vector of initial guesses for the parameters
%Setup to fit the data with the solution of the ODE
options=optimset('FunValCheck','on','MaxIter',5000,'MaxFunEvals',100000,'TolX',1e-8,'Tolfun',1e-8);
[pfinal,resnorm,exitflag,output]=lsqcurvefit(@fitN,p,xdata,ydata,[],[],options);
%Plot output of lsqcurvefit
xsim = linspace(xmin,xmax,50);
yout = fitN(pfinal,xsim);
exitflag
output
pfinal
%plot the solution to the ODE from bvp4c
figure(1)
hold on
plot(xdata,ydata,'go');
xlabel('Depth')
ylabel('Nitrate (µol/cm3)')
legend('Nitrate model','Nitrate data')
%connects bvp4c to lsqcurvefit
plot(xdata,ydata,'go');
xlabel('Depth')
ylabel('Nitrate (µol/cm3)')
legend('Nitrate model','Nitrate data')
function F = fitN(x,xdata)
global Gm N0 xhat p xmin xmax Dno3
p=x;
options= bvpset('Stats','on','RelTol',1e-4,'abstol',1e-3);
solinit = bvpinit(linspace(xmin,xmax,50),[0.022,0.03]);
sol = bvp4c(@ode, @bc, solinit,options);
ydata_sim = deval(sol,xdata)'; %evaluates solution of bvp at the xdata points
F = ydata_sim(:,1);
%The ODE for Nitrogen
function dydx = ode(x,y)
global Gm N0 xhat p xmin xmax Dno3
f=zeros(length(x),1);
for i=1:length(x)
if (x(i) <= xhat)
A = 1;
else
A = 0;
end
end
%The boundary conditions for Nitrogen
function boundary = bc(ya,yb)
global Gm N0 xhat p xmin xmax Dno3
Torsten
Mar 23, 2012, 3:56:16 AM3/23/12
to Lauren
1. I'm not sure about MATLAB conventions, but at least in FORTRAN the
ordering and number of variables in global statements must be the same
throughout the complete coding. So for security reasons you may delete
the variable "g" from the first global statement.
2. To account for the variable "Dno3", you should make the fitting without
this variable (thus with the "old" differential equation) and multiply
parameters p1, p3 and p5 afterwards by Dno3. This is equivalent to
fitting the "new" differential equation.
Best wishes
Torsten.
0 new messages