How do you make an initial guess in bvp4c?

2 Ansichten (letzte 30 Tage)
Marcus Rosales
Marcus Rosales am 5 Feb. 2024
Kommentiert: Marcus Rosales am 6 Feb. 2024
I am trying to sovle the following differential equation using bvp4c:
There are issues at the origin... I need to solve here though!
In certain limits you can say , so I would like to guess this as an initial solution. I tried something like this (xl=0, xr=100 or so):
function [Sxint, xint] =bvp4c_mathworks_Abrikosov_guess(xl,xr)
N = 5000;
xint = linspace(xl,xr,N);
solinit = bvpinit(xint,@guess);
sol4c = bvp4c(@myode,@mybc,solinit);
Sxint = deval(sol4c,xint);
end
function dy = myode(r,y)
k = 5/sqrt(2);
dy(1,1) = y(2);
dy(2,1) = k^2*(y(1)^2-1)*y(1)+y(1)/r^2-y(2)/r;
end
function res = mybc(ya,yb)
res = [ya(1)
yb(1)-1];
end
function g = guess(x)
k = 5/sqrt(2);
g = tanh(x/k);
end
Is there a way to do this? What do I need to change?

Akzeptierte Antwort

Torsten
Torsten am 5 Feb. 2024
Bearbeitet: Torsten am 5 Feb. 2024
xl = 1e-8;
xr = 5;
[Sxint, xint]=bvp4c_mathworks_Abrikosov_guess(xl,xr);
plot(xint,Sxint(1,:))
hold on
M = cell2mat(arrayfun(@(x)guess(x),xint,'UniformOutput',0));
plot(xint,M(1,:))
hold off
function [Sxint, xint] =bvp4c_mathworks_Abrikosov_guess(xl,xr)
N = 500;
xint = linspace(xl,xr,N);
solinit = bvpinit(xint,@guess);
sol4c = bvp4c(@myode,@mybc,solinit);
Sxint = deval(sol4c,xint);
end
function dy = myode(r,y)
k = 5/sqrt(2);
dy(1,1) = y(2);
dy(2,1) = k^2*(y(1)^2-1)*y(1)+y(1)/r^2-y(2)/r;
end
function res = mybc(ya,yb)
res = [ya(1)
yb(1)-1];
end
function g = guess(x)
k = 5/sqrt(2);
g = [tanh(x/k);1/k*(1-tanh(x/k)^2)];
end
  3 Kommentare
Torsten
Torsten am 5 Feb. 2024
Yes, you solve for f and f', and you have to specify initial profiles for both of them.
So if you specify f0 = tanh(x/k) as initial guess for f, it's natural to use f0' = 1/k*(1-tanh(x/k)^2) as initial guess for f'.
Marcus Rosales
Marcus Rosales am 6 Feb. 2024
Makes sense, it needs a way of knowing how the function changes under a step. Thanks again!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte


Version

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by