Solve BVP with datapoints as input variable

Hello,

I have been trying to solve first order coupled differential equations in which both of my differential equations have variables (alpha_e and alpha_h) as an input. The variables are simulated results, and thus I have values not an equation to represent these variables. When I generate a fit equation and use them directly then I can solve it. But finding an equation to fit the data is not always possible, so I need to use the values directly/some other way (defining a function using fit command), but nothing seems to work except for equation. I am getting 'Not enough input arguments' when I am trying with fit function. I am using BVP5c solver. Please help.

Here are the equations and the codes that I am working with. If you run the code, you will see the expected results, but if you use the ones defined by 'fit' (bold) will return: 'Not enough input arguments'

Codes:

solinit = bvpinit(y, @guess);

sol = bvp5c(@odefun, @bcfun, solinit);

function dpdy = odefun(y,p, e,h, ealpha, halpha)

% ealpha = fit(y, e, 'linearinterp', 'Normalize', 'on');

% halpha = fit(y, h, 'linearinterp', 'Normalize', 'on');

ealpha = ((140800*exp(-((y-0.17)/0.008511).^2)+112300*exp(-((y-0.1832)/0.01434).^2)-186.6*exp(-((y-0.06409)/1.297).^2)+3531000*exp(-((y-0.2225)/0.02767).^2)+70120*exp(-((y-0.2767)/0.05906).^2)-3416000*exp(-((y-0.2227)/0.02722).^2)));

halpha = ((53050*exp(-((y-0.1616)/0.009634).^2)-12310*exp(-((y+42.82)/18.3).^2)+12320*exp(-((y-0.2441)/0.02159).^2)+49680*exp(-((y-0.1774)/0.01608).^2)+50640*exp(-((y-0.2046)/0.02749).^2)+33.49*exp(-((y-0.06427)/0.04413).^2)+21050*exp(-((y-0.2638)/0.05081).^2)));

dpdy = zeros(2,1);

dpdy(1) = (1-p(1)).*(ealpha./10000).*(p(1)+p(2)-p(1).*p(2));

dpdy(2) = -(1-p(2)).*(halpha./10000).*(p(1)+p(2)-p(1).*p(2));

end

function res = bcfun(pa,pb)

res = [pa(1)

pb(2)];

end

function p = guess(y) % guess at solution behavior

p = [sin(y)

cos(y)];

end

[Torsten]

You don't need to use bvp5c.

The solution is obviously P_e = P_h = 0 for all y.

[Mohammad Rezaul Islam]

No. This is the solution. Why do you think that?

[Torsten]

No. This is the solution. Why do you think that?

Because P_e = P_h = 0 solves your differential equations and your boundary conditions.

Maybe there is more than one solution, but usually in such cases there are infinitly many, and numerical methods are senseless to apply.

[Mohammad Rezaul Islam]

Thank you, Torsten for your replies.

You can see the solution for yourself, just add this line to the code: plot(sol.x,sol.y);

There is one unique solution depending on the values of alpha_e and alpha_h, as far as I understand.

However, from Matlab/coding point of view, can anyone help me with a proper definition of ealpha and halpha that will allow the 'odefun' to execute without an error and producing a result?

[Torsten]

There is one unique solution depending on the values of alpha_e and alpha_h, as far as I understand.

Then you didn't understand that P_e = P_h = 0 is always a solution - whatever ealpha and halpha are. So there is no unique solution for your system.

However, from Matlab/coding point of view, can anyone help me with a proper definition of ealpha and halpha that will allow the 'odefun' to execute without an error and producing a result?

Use MATLAB's "interp1" with method = spline on your data for alpha_e and alpha_h instead of "fit".

[Mohammad Rezaul Islam]

Use MATLAB's "interp1" with method = spline on your data for alpha_e and alpha_h instead of "fit".

I have tried, same error.

"Not enough input arguments.

Error in BVP_PePh_usingInterpolation>odefun (line 16)

ealpha = interp1(y, e, 'spline');

Error in bvparguments (line 105)

testODE = ode(x1,y1,odeExtras{:});

Error in bvp5c (line 135)

bvparguments(solver_name,ode,bc,solinit,options);

Error in BVP_PePh_usingInterpolation (line 11)

sol = bvp5c(@odefun, @bcfun, solinit);"

[Mohammad Rezaul Islam]

Dear Torsten,

Thanks a lot for solving the problems in my code, your responses were very fast and effective. The code runs without an error and produces results.

However, using 'interp1' listed methods, the result varies and I am not getting an expected result. I did check the data for errors, and removed unnecessary inputs (values 1E-34.. etc to 0), but still gives different result than when I am using Fit Equations. If you have any solution/insight to this, please share, I would be very happy if you can help me further.

Additionally:

Because P_e = P_h = 0 solves your differential equations and your boundary conditions.

This is a trivial solution, we are looking for a non-zero result for a certain range of input with certain coefficient values, the equation is not valid for all values of y, as far as I know, according to the theory.

[Mohammad Rezaul Islam]

Might be the problem is from the simulated data. I am working on this.

Thanks.

[Torsten]

One problem might be that the input data are not smooth enough, but since they are simulated data, they usually should work for simply interpolating them between the y-values.

[Mohammad Rezaul Islam]

Yes. The values from my recent simulations include one of the part of the ODEs into account ((p(1)+p(2)-p(1).*p(2)), as a results, the simulated ealpha_array and halpha_array are different than before. By comparing the new data with the old one, I have found that. If I move back to my earlier simulation settings, then I am able to get expected result.

I am so grateful to you for helping me out. Wish you all the best, Torsten. The problem is successfully solve with your essential help.

Reza

Answer 1

[Torsten]

Pass the arrays of y-values, e_alpha and h_alpha to odefun:

y_given = some 1d-array;
ealpha_given = some 1d-array of the same size as y_given;
halpha_given = some 1d-array of the same size as y_given;
sol = bvp5c(@(y,p)odefun(y,p,y_given,ealpha_given,halpha_given), @bcfun, solinit);

and interpolate to the y-values from bvp5c:

function dpdy = odefun(y,p, y_given,ealpha_given, halpha_given)
  ealpha = interp1(y_given,ealpha_given,y,'spline');
  halpha = interp1(y_given,halpha_given,y,'spline');
  dpdy = zeros(2,1);
  dpdy(1) = (1-p(1))*(ealpha/10000)*(p(1)+p(2)-p(1)*p(2));
  dpdy(2) = -(1-p(2))*(halpha/10000)*(p(1)+p(2)-p(1)*p(2));  
end