Using a given array as an initial guess for bvp4c using bvpinit
7 views (last 30 days)
Show older comments
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
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.
Answers (0)
See Also
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!