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

6 次查看（过去 30 天）

Marcus Rosales 2024-5-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2118836-using-a-given-array-as-an-initial-guess-for-bvp4c-using-bvpinit

评论： Torsten 2024-5-15

bvp4c_mathworks_half_VAV_2.m

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 个评论
显示 1更早的评论隐藏 1更早的评论

Marcus Rosales 2024-5-15

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).