Using a given array as an initial guess for bvp4c using bvpinit

7 views (last 30 days)
I tried using the solution, G, to the "free" version of a given ODE (given in code with a=0) to initialize the bvp4c solver. Given a mesh xint and a function guess(x) (made by interp1, G and xint), I used the following:
solinit = bvpinit(xint, @guess);
The details are in the attached code (forgot to include: xl=-50, xr=2). The free case works with @guess replaced with [0 0].
I keep getting the error: "Index exceeds matrix dimensions" when the solver attempts the ode solution. I have tried a lot to get this to work, but things keep failing. Any clues on how I am probbaly screwing up "solinit"?
function [Sxint, x] =bvp4c_mathworks_half_VAV_2(xl,xr, N, theta, r0, Guess)
global tt
global r
global G
global xx
xint = linspace(xl,xr,N); %BC at [-inf, inf] replaced with [xl, xr]
xx = xint;
tt = theta; %parameters: winding number and vortex radius, say pi/6, 16
r = r0;
G = zeros(1,length(Guess)); %What I want to initialize with
solinit = bvpinit(xint, @guess);
sol = bvp4c(@myode,@mybc,solinit);
Sxint = deval(sol,xint);
end
function dy = myode(t,y)
size(y)
size(t)
a = 0*10^(-5); %a reduced coupling to pdw times uniform gap here
dy(1,1) = y(2);
dy(2,1) = 1/2*(1-exp(2*t)*cos(2*y(1)))*sin(2*y(1))+a*a*exp(2*t)*cos(f(t))*cos(2*y(1));
end
function phase_diff = f(t)
global tt;
global r;
dd = 80;
phase_diff = pi*tanh(r*exp(t)/dd)-tt;
end
function res = mybc(ya,yb)
res = [ya(1)
yb(1)-pi/4];
end
function [y1, y2] = guess(x)
global G;
global xx;
ee = 10^(-10);
y1 = interp1(xx,G,x).';
y2 = (interp1(xx,G,x+ee).'-y1)/ee;
%y2 = gradient(G).';
end
  3 Comments
Marcus Rosales
Marcus Rosales on 15 May 2024
I originally used the gradient here. Will this also give me some issues? I figured the end point might run into some issues since there is nothing to interpolate with (should be zero though).
Torsten
Torsten on 15 May 2024
If you can supply G, don't you have access to dG/dx ?
But the bad approximation for dG/dx is not the reason for the error message. As said: If xx and G have the same length, your code should work.

Sign in to comment.

Answers (0)

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!