Solve DET and ODE faster

[EDIT: 20110617 14:03 CDT - reformat - WDR]

Hye I made a set of function, but it is too slow or innacurate. Here are the functions:

first one-------------------------------------------

function Y=plasmainhr(s,q,a,b,Ed,f,mu1,u,N,x0)
F=Deter(s(1),q,a,b,Ed,f);
x=fsolve(F,x0);
for j=1:length(s)
    if j==1
        Y(j)=EDOplasma(s(j),q,a,b,Ed,f,x,mu1,u,N);
    else
        Y(j)=EDOplasma(s(j),q,a,b,Ed,f,Y(j-1),mu1,u,N);
    end
end

for k=1:length(Y)
    beta(k)=real(Y(k));
    alpha(k)=abs(imag(Y(k)));
end

subplot(2,1,1)
plot(beta,s,'k')
xlabel('beta(1/m)')
ylabel('omega/omegap')
title('BETA')
subplot(2,1,2)
plot(alpha,s,'k')
xlabel('alpha(1/m)')
ylabel('omega/omegap')
title('ALPHA')

2nd------------------------------------------------------

function X=EDOplasma(s,q,a,b,Ed,f,x0,mu1,u,N)
beta=zeros(1,N+2);
beta(2)=x0;
omega=f*2*pi;
v=q*omega;

Eo=8.85418782*10^-12;
muo=1.25663706*10^-6;
me=9.10938188*10^(-31);
e=-1.602176487*10^(-19);

ne0=omega^2*me*Eo*a^2/(2*s^2*e^2*(quadl(@(x)(besselj(0,mu1.*x./a).*x),0,a)));

Ko=sqrt((omega^2)*Eo*muo);

n=1;
beta(1)=0;
while abs(beta(n)-beta(n+1))>u
    n=n+1;
    if n-2==N
        fprintf('La convergence désiré de:%-3.2e n''a pas pu être atteinte en %i itérations',u,N)
        beta(n+1)=beta(n);
        break
    end
    c=beta(n);
    save('B','c','ne0','Ko','a','mu1','omega','v')
    opts = odeset('AbsTol',5e-2);
    [t,y]=ode23(@Gi,[0.000001;a],[1,0],opts);
    F=@(x)10000000000000*det([0 besselj(0,(sqrt(Ko^2*Ed-(x)^2))*b) (i*bessely(0,sqrt((Ko^2*Ed-(x)^2))*b)+besselj(0,sqrt((Ko^2*Ed-(x)^2))*b)) -(i*bessely(0,sqrt((Ko^2-(x)^2))*b)+besselj(0,sqrt((Ko^2-(x)^2))*b));y(length(t),1) -besselj(0,(sqrt(Ko^2*Ed-(x)^2))*a) -(i*bessely(0,sqrt((Ko^2*Ed-(x)^2))*a)+besselj(0,sqrt((Ko^2*Ed-(x)^2))*a)) 0;0 -Ed*besselj(1,(sqrt(Ko^2*Ed-(x)^2))*b)/((sqrt(Ko^2*Ed-(x)^2))) -Ed*(i*bessely(1,sqrt((Ko^2*Ed-(x)^2))*b)+besselj(1,sqrt((Ko^2*Ed-(x)^2))*b))/((sqrt(Ko^2*Ed-(x)^2))) (i*bessely(1,sqrt((Ko^2-(x)^2))*b)+besselj(1,sqrt((Ko^2-(x)^2))*b))/((sqrt(Ko^2-(x)^2)));-(1-((e^2*ne0*besselj(0,mu1)/(me*Eo))/(omega*(omega-i*v))))*-y(length(t),2)/(((Ko^2*(1-((e^2*ne0*besselj(0,mu1)/(me*Eo))/(omega*(omega-i*v))))-(x)^2))) Ed*besselj(1,(sqrt(Ko^2*Ed-(x)^2))*a)/((sqrt(Ko^2*Ed-(x)^2))) Ed*(i*bessely(1,sqrt((Ko^2*Ed-(x)^2))*a)+besselj(1,sqrt((Ko^2*Ed-(x)^2))*a))/((sqrt(Ko^2*Ed-(x)^2))) 0]);
    beta(n+1)=fsolve(F,beta(n));
end
X=beta(n+1);

3rd------------------------------------------

function rk=Gi(t,y)

Eo=8.85418782*10^-12;
muo=1.25663706*10^-6;
me=9.10938188*10^(-31);
e=-1.602176487*10^(-19);
load B;

rk(1,1)=y(2);
rk(2,1)=(-((1/t)-c^2*e^2*ne0*mu1*besselj(1,mu1*t/a)/((Ko^2*(1-(e^2*ne0*besselj(0,mu1*t/a)/(me*Eo*omega*(omega-i*v))))-c^2)*(1-(e^2*ne0*besselj(0,mu1*t/a)/(me*Eo*omega*(omega-i*v))))*me*Eo*omega*(omega-i*v)*a))*y(2)-(Ko^2*(1-(e^2*ne0*besselj(0,mu1*t/a)/(me*Eo*omega*(omega-i*v))))-c^2)*y(1));

4th-----------------------------------------------

function F=Deter(s,t,a,b,Ed,f)

omega=f*2*pi;

v=t*omega;

omegap=omega/s;

Eo=8.85418782*10^-12;
muo=1.25663706*10^-6;

Ko=sqrt((omega^2)*Eo*muo);
Ep=1-((omegap^2)/(omega*(omega-i*v)));

x=sym('x');

F=@(x)1000000*det([0 besselj(0,(sqrt(Ko^2*Ed-(x)^2))*b) (i*bessely(0,sqrt((Ko^2*Ed-(x)^2))*b)+besselj(0,sqrt((Ko^2*Ed-(x)^2))*b)) -(i*bessely(0,sqrt((Ko^2-(x)^2))*b)+besselj(0,sqrt((Ko^2-(x)^2))*b));besselj(0,(sqrt(Ko^2*Ep-(x)^2))*a) -besselj(0,(sqrt(Ko^2*Ed-(x)^2))*a) -(i*bessely(0,sqrt((Ko^2*Ed-(x)^2))*a)+besselj(0,sqrt((Ko^2*Ed-(x)^2))*a)) 0;0 -Ed*besselj(1,(sqrt(Ko^2*Ed-(x)^2))*b)/((sqrt(Ko^2*Ed-(x)^2))) -Ed*(i*bessely(1,sqrt((Ko^2*Ed-(x)^2))*b)+besselj(1,sqrt((Ko^2*Ed-(x)^2))*b))/((sqrt(Ko^2*Ed-(x)^2))) (i*bessely(1,sqrt((Ko^2-(x)^2))*b)+besselj(1,sqrt((Ko^2-(x)^2))*b))/((sqrt(Ko^2-(x)^2)));-Ep*besselj(1,(sqrt(Ko^2*Ep-(x)^2))*a)/((sqrt(Ko^2*Ep-(x)^2))) Ed*besselj(1,(sqrt(Ko^2*Ed-(x)^2))*a)/((sqrt(Ko^2*Ed-(x)^2))) Ed*(i*bessely(1,sqrt((Ko^2*Ed-(x)^2))*a)+besselj(1,sqrt((Ko^2*Ed-(x)^2))*a))/((sqrt(Ko^2*Ed-(x)^2))) 0]);

--------------------------------------

To call them we use the line:

Y=plasmainhr(s,q,a,b,Ed,f,mu1,u,N,x0)

q=0.001; a=0.013; b=0.015; Ed=4.52; f=200000000; mu1=0; N is the maximum number of iteration in the second function; u is the accuracy we ask to end the second function and x0=8.5 for s(1)=0.1; s is the vector we use to plot and find the solutions. I usually use 0.1:0.0001:0.18, but I should up to 0.4.

The problem is that the second function take too much time. I changed u, N, the increment in s, the accuracy in odeset, but when I go fast enough the anwers diverge. The graph should be a smooth curved line.

Thanks.

[Walter Roberson]

Please check your second function; the second last line is obviously only a partial statement, as if something was not pasted in correctly.

[Liber-T]

I've corrected it.

Answer 1

[Walter Roberson]

Is there no way to make a guess at the starting point, and no way to put reasonable bounds on the search interval? fsolve() has much better performance if you can put good bounds on the search.

Sometimes calculations are just slow.

It appears to me that there are some expressive sub-expressions in your F function that occur several times (e.g., multiplied by different things.) You might be better off not using an anonymous function and instead coding a more optimized version of the function.

Is this det() the same long function you were looking at before, that Matt Fig and I worked with you on a bit? If it is I could go back to that one and run it through an optimizer.

[Liber-T]

It appears to me that there are some expressive sub-expressions in your F function that occur several times (e.g., multiplied by different things.) You might be better off not using an anonymous function and instead coding a more optimized version of the function.

[Walter Roberson]

I bet you tried to quote what I said using "" around what I wrote... Unfortunately there is a bug in MATLAB Answers that causes the rest of your response to disappear when you do that! If this happened to you, please go to http://www.mathworks.com/matlabcentral/answers/994-wish-list-for-matlab-answer-sections

and vote on either my Answer that starts "When constructing a comment", or on Matt Fig's Answer that starts "I am adding this again because it is REALLY, VERY annoying". (Or both!)

[Liber-T]

So, a summary of what I said:

How do I bound fslove?

I'm not sure to understand what you mean by the thing I quote and erased the rest of my text.

And the determinant here is a little bit different from the other one, and the det in the 2nd and 4th function are different from each other. If you want to rearrenge them it's alrigth as long as mathmatically they stay the same.

Thanks

[Walter Roberson]

I keep thinking that fsolve() is like Maple's fsolve(), but it isn't. There does not appear to be any mechanism to specify bounds for MATLAB's fsolve().

Possibly you might do better with a different 'Algorithm' option; scroll down to the "Options" section of the fsolve() documentation.

If I recall correctly, you have the symbolic toolbox; if so then you could try using numeric::solve() with the RestrictedSearch option; see http://www.mathworks.com/help/toolbox/mupad/numeric/solve.html

Wait... you are searching for a real root of an continuous function of one variable, aren't you? If so then you could try using fzero() instead of fsolve(). fzero() allows a search interval to be specified.

[Liber-T]

In fact it's a complex root of a complex function. Normally used, fzero and solve doesn't work.