ode45/fsolve issue for Shooting Method solution to boundary value problem

I'm trying to solve a boundary value problem in MATLAB using the shooting method. However, when I try to pass my function through fsolve, I'm getting errors like: Warning: Failure at t=-9.462647e+001. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.273737e-013) at time t.
I've tried some of the common workarounds like using different ODE solvers and changing the tolerances in the ODE solve and fsolve but to no avail. Currently, when I run fsolve on my parameters, it changes nothing and it's the same as my initial guess.
I've included the code for the main script file as well as the function file called by fsolve.
%%Main file
% PARAMETERS
% propanol (fig2, curve3, p337)
E=6.5; %2.5
K=0.0474;
% GUESSES
infinity=100; % numerical approximation of infinity
x=0; % expanded around x=0
x0=3;
B0=-2*10^(-8);
B1=0.01;
% Guess for P
P=-1+3*B0*(1-(E/(K+1)))*exp((x-x0)*sqrt(3*E/(K+1)));
% Guess for P'
PI=sqrt(3*E/(K+1))*3*B0*(1-(E/(K+1)))*exp((x-x0)*sqrt(3*E/(K+1)));
% Guess for H
H=1+B0*exp((x-x0)*sqrt(3*E/(K+1)))+B1*exp((x-x0)*sqrt(3));
% Guess for H' (contact angle in degrees)
HI=sqrt(3*E/(K+1))*B0*exp((x-x0)*sqrt(3*E/(K+1)))+sqrt(3)*B1*exp((x-x0)*sqrt(3));
% FUNCTION
f=@(t,x) [x(2); ...
3*E*(1+x(1))/((1+K*x(3))*(x(3))^3) - 3*x(2) * x(4)/(x(3));...
x(4);...
-1/(x(3)^3) - x(1)];
% SHOOTING METHOD
tspan=[-infinity, infinity];
myguess=[P, PI, H, HI]; % initial guess for P', H'
[actualvalue, Fval, exitflag]=fsolve('shootingfunc_thin_100',myguess);
init=actualvalue; % new initial paramaters
[t,X]=ode45(f,tspan,init);
%%Function file
function [ f ] = shootingfunc_thin_100( a )
%shootingfunc_thin: Shooting from thin to thick
% optimizes P=0 at infinity only
E=6.5;
K=0.0474;
fun=@(t,x) [x(2); ...
3*E*(1+x(1))/((1+K*x(3))*(x(3))^3) - 3*x(2) * x(4)/(x(3));...
x(4);...
-1/(x(3)^3) - x(1)];
infinity=10^2;
tspan=[-infinity, infinity];
init=[-1 a(2) 1 a(4)]; % initial conditions at thin
%options=odeset('RelTol',1e-10); %'RelTol',1e-16
[t,x]=ode45(fun,tspan,init);
f=x(end,1);
end

回答 (0 件)

カテゴリ

質問済み:

2012 年 11 月 28 日

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by