OCTAVE - Equazione trascendente

[guido.bonalumi]

Ciao a tutti,
innanzitutto mi scuso nel caso avessi postato nella sezione sbagliata; se così fosse vi pregherei di indirizzarmi a quella giusta.
Venendo al mio problema... sto risolvendo un piccolo problema di scambio termico la cui soluzione richiede di ricavare le prima N (devo ancora capire bene quante) radici della seguente equazione:

x*(J1(x))/(J0(x))=Bi

dove Bi è il numero di Biot (Bi=h*L/k, che possiamo assumere pari a 0,02) e J0 e J1 sono due funzioni di Bessel del primo tipo (Jn(n, x)). Sono in grado di risolvere le funzioni di Bessel separatamente (besselj(n, x)), ma non ho la minima idea di come fare per estrarre le radici di cui ho bisogno dall'equazione che ho scritto sopra.
Qualcuno saprebbe come potrei fare?

Vi ringrazio anticipatamente per l'aiuto che saprete darmi.

[fibonacci]

puoi usare il calcolo simbolico in Matlab

[guido.bonalumi]

Grazie per la risposta, ma potresti essere un pò più dettagliato? Non sono molto esperto ne di octave ne tantomeno di matlab (anche se so che sono praticamente uguali).

[fibonacci]

allora non ne sono certo al 100 per cento ma non si può fare calcolo simbolico in octave
per questo ti ho detto di usare MATLAB
è davvero semplice da usare
esempio
syms x a
eq='sin(x)+a=0'
solve(eq,x)
ti verranno restituite le soluzioni in funzione del parametro a
cioè in solve hai indicato di risolvere l'equazione rispetto alla variabile x

[guido.bonalumi]

Ti ringrazio. Il problema è che non ho a disposizione MATLAB per questo uso octave. Ho provato con la funzione fsolve in questo modo:

Codice: Seleziona tutto

fsolve("(besselj(1, x))/(besselj(0, x))-Bi/x", 1.1)

ma mi restituisce il seguente errore:

Codice: Seleziona tutto

error: @(besselj(1, x))/(besselj(0, x))-Bi/x: no function and no method found
error: called from:
error:   /usr/share/octave/3.6.1/m/optimization/fsolve.m at line 150, column 9
error:   /home/guido/Documenti/Università/trasmissione del calore/progetto/octave/transcendental_eqn.m at line 7, column 4

Non so se questa funzione dovrebbe fare qualcosa di simile a quella che mi hai indicato tu.

[fibonacci]

prova con questa sintassi
come lo hai usato tu sicuro non funziona perché non è un programma di calcolo simbolico Octave

f = @(x) (0.2)^x-log(x)
fsolve(f,1.1)

naturalmente devi mettere la tua funzione questo solo un esempio per la sintassi

[bite]

Questa equazione con le funzioni di Bessel non puoi risolverla simbolicamente, devi andare con il metodo di Newton o qualcosa del genere.

Puoi usare Wolfram Alpha. Vai su http://www.wolframalpha.com, schiaffa la stringa

Codice: Seleziona tutto

x BesselJ[1, x] - 0.02 BesselJ[0, x] == 0

nella casellina di testo e clicca il segno di uguale.

[guido.bonalumi]

Vi ringrazio per le risposte. Con la soluzione di fibonacci funziona! Ho risolto in questo modo:

Codice: Seleziona tutto

f = @(x) (x*(besselj(1, x))/(besselj(0, x))-Bi)
fsolve(f,1.1)

Solo un'ultima cosa... io avrei bisogno delle prime N radici, mentre in questo modo ne trovo solo una, quella prossima a 1,1. Come potrei fare per evitare di doverli cercare manualmente? Stavo pensando un'iterazione ma se non conosco i valori per cui approssimare non ho idea di come fare. Qualche idea?

[guido.bonalumi]

Ho quasi risolto... sono giunto a questo:

Codice: Seleziona tutto

D=0.05
h=100
k=60.5
Lc=D/4
Bi=h*Lc/k
N=10;

# Plot della funzione nell'intervallo compreso fra x=0 e x=20 per avere un'idea dei valori delle radici ricercate
x=[0.0001:0.001:20];
y=x.*(besselj(1, x))./(besselj(0, x)).-Bi;
plot(x, y, '-r')
axis([0 20 -1 1])

a=[1:N];
for i=0:(N-1)
	f=@(x)(x*(besselj(1, x))/(besselj(0, x))-Bi);
	a(i+1)=fsolve(f, (i*pi+pi/4));
end
a

Il risultato di per a è il seguente:

Codice: Seleziona tutto

a =

    0.20276    3.83709    7.01853   10.17550   13.32524   16.47188   19.61691   22.76099   25.90447   29.04754

Il che è giusto solo parzialmente. Infatti se si prova a guardare il grafico della funzione, si può notare come i valori salvati nel vettore a siano solo le soluzioni dispari. Fra ogni coppia di valori salvati infatti ve n'è un altro che non viene salvato.
Ho provato a raffinare al soluzione passando da

Codice: Seleziona tutto

i*pi+pi/4

a

Codice: Seleziona tutto

i*pi/2+pi/4

ma il risultato è il seguente:

Codice: Seleziona tutto

a =

    0.20276    0.20276    3.83709    3.83710    7.01853    7.01853   10.17550   10.17551   13.32524   13.32526

Evidentemente non è quello che voglio, infatti vengono trovate le stesse radici di prima, semplicemente ciscuna viene trovata due volte. C'è qualcosa che mi sfugge, probabilmente riguardo al funzionamento di fsolve. Qualcuno saprebbe dirmi cosa sbaglio?

[bite]

Il grafico della funzione x BesselJ[1,x]/BesselJ[0,x]-0.02 è questo:



Tra le radici che hai elencato non ce ne sono altre, ci sono solo dei punti di discontinuità di seconda specie in corrispondenza degli zeri di BesselJ[0,x].

[BlueEyes]

Confermo, l'algoritmo "salta" metà delle radici, quelle indicate nel grafico con le frecce blu. Ciao

[bite]

Ripeto: quelle non sono radici ma punti di discontinuità di seconda specie.

Che le radici siano solo quelle elencate lo si può vedere plottando la funzione x BesselJ[1, x] - 0.02 BesselJ[0, x] che non ha punti di discontinuità ed ha le stesse radici:

[BlueEyes]

Retromarcia, la mia!
Si vede chiaramente, modificando leggermente la funzione, che quelli intermedi sono asintoti, o comunque segni grafici di gnuplot quando la funzione attraversa le discontinuità. Come da grafico e codice allegati. Ciao

Codice: Seleziona tutto

N=3;
x=[0:0.1:5];
y=(besselj(1, x))./(besselj(0, x));
plot(x, y, '-r')
axis([0 5 -1 1])
grid

[bite]

Sono asintoti verticali