Using a given array as an initial guess for bvp4c using bvpinit
7 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Torsten
el 15 de Mayo de 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.
Respuestas (0)
Ver también
Categorías
Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!