Equazione differenziale con metodo di Newton

[feddy]

Sbagli dal principio, quella è una ODE, hai solo i dati a $x=0$, non hai nessuna condizione sull'estremo destro. Non puoi discretizzare l'intervallo se un intervallo non ce l'hai. Questa è roba base base, ti conviene prenderti (scaricare) un libro sulla risoluzione numerica di equazioni differenziali. O magari qualche dispensa tipo questa, ma è solo un esempio.
Devi ricondurti a un sistema del prim'ordine introducendo delle variabili ausiliarie e poi applicare uno schema numerico opportuno (Eulero implicito, trapezi, Runge Kutta). Se non vuoi scriverti il codice del metodo basta prendere MatLab e usare odexyz e trovarti la soluzione.

Quello che segue era la mia prima risposta a questo problema. Non mi ero accorto che fosse una ODE e mi sono concentrato solo sul risponderti sulla condizione di Robin del BVP.

EDIT

Tutto quello scritto sotto vale per la condizione di Robin nel caso di problemi ai limiti, ma non risolve l'ultimo problema dell'OP.

Ora per me sono giorni un po' pesanti (come lo sono per te immagino), ma appena ho un po' di tempo do' meglio una occhiata al tuo codice.

Per la condizione di Robin a $x=0$ provo a dirti come farei io, senza passare dallo sviluppo di Taylor, per comodità:

Questa è $\frac{d}{dx}u(0) + \frac{q_0}{k_0 + k_1 u+k_2 u^2}=0$

Discretizzo questa riga
$\frac{u_0 -2u_1 +u_2}{h^2} + \frac{g}{k} =0$, ($\square$) dove con $u_0$ indico il valore della soluzione nel nodo fantasma, mentre $h$ è il passo spaziale. Chiaramente il problema ora è la ricerca del valore della soluzione al nodo fantasma, e questo viene ricavato dalla condizione di Robin stessa, che ora scrivo per esteso:

$\frac{u_2 - u_0}{2h} + \frac{q_0}{k_0 + k_1 u_1 + k_2 * u_{1}^{2}}=0$

Ora da qui ricavi $u_0$, lo sostituisci in $\square$ e ricavi l'espressione della prima riga, che poi modificherai di conseguenza

Fine EDIT

[luca.milano]

Ma figurati, sei stato gentilissimo, prenditi pure il tempo che ti serve per fare i tuoi esami e più avanti se avrai voglia potrai ripassare per questa discussione. In ogni caso mi hai già aiutato tanto e per questo ti ringrazio!

Purtroppo il tuo ultimo commento ammetto di non averlo capito benissimo. In particolare, lo stesso problema ma svolto in modo lineare, quindi con $ k $ costante nota, riesco a farlo funzionare e forse per questo pensavo fosse relativamente semplice estenderlo al caso non lineare utilizzando il metodo di Newton.

Giusto per completezza ti riporto il codice del caso lineare, che magari guarderai più avanti quando avrai più tempo:

Codice:

clear
clc
%CODICE CAMPO DI T CON ASSEGNATA TEMPERATURA E FLUSSO IN INGRESSO

L=10;
N=100; %numero elementi
dx=L/N;

%SPECIFICA LA TEMPERATURA ALL'INGRESSO E IL FLUSSO TERMICO
T0=100;
q0=-500;

%COEFFICIENTI DELL'EQUAZIONE: generazione di potenza e conduttività termica
g=250;
k=50; %conduttività termica

A=zeros(N,N);
b=zeros(N,1);

A(1,1)=1; %prima riga
A(2,2)=+1; %seconda riga

b=[T0 (T0-dx/k*q0-1/2*dx^2*g/k)]';
for j=3:N
A(j,j-2)=1;
A(j,j-1)=-2;
A(j,j)=1;

b(j,1)=-g/k*dx^2;
end

T=A\b;

x=linspace(0,L,N);
plot(x,T)


Lo avevo scritto prima che tu mi spiegassi come scrivere le matrici agilmente, quindi ci sono ancora i cicli for.

Buona fortuna per i tuoi esami!

[feddy]

Buona fortuna per i tuoi esami!

Grazie !

Ad ogni modo ti ripeto che per risolvere ODE, che è quello che vuoi fare, devi ricondurre l'equazione diff. (che era al secondo ordine) in un sistema al prim'ordine e usare MatLab. Il fatto che plottasse qualcosa non vuol dire niente in realtà. Ad esempio, tu hai messo esplicitamente una condizione sul bordo a $x=L$. Quindi, faccio finta che sulle ascisse ci sia il tempo e ti scrivo 5 righe di come andrebbe fatto uno stupidissimo Eulero Esplicito per esempio col tuo problema

\( \begin{cases} (d^2 u)/dx^2 +g/k=0 \\ u(0)=u_0\\ \frac{du}{dx}(0)=-\frac{q_0}{(k(u_0))} \end{cases} \)

Ho scritto il problema al prim'ordine introducendo $u_{1}'=u_2$ e risolto un sistema differenziale con condizioni iniziali che ho scelto io, visto che tu precedentemente non le davi (ho solo assunto $u(0)=u_0=0$ poiché suppongo che la temperatura all'istante iniziale sia nulla, comunque non cambia niente, basta che sistemi $y_0$). Ancora, visto che il risultato comparirà sulla tua tesi, se vuoi garantirti un risultato migliore, senza preoccuparti della tanta teoria che ci sta dietro, puoi usare il comando ode di MatLab, che mi pare tu conosca già.

Ad ogni modo, i tuoi problemi sono più che altro sulla teoria, e continuare a scrivere pezzi di codice di problemi diversi rende questa discussione veramente "enorme" e si rischia di perdere di vista gli obiettivi. Ti ripeto che secondo me ti conviene darti una letta a del materiale sull'argomento, almeno saper distinguere una ODE da un problema ai limiti, altrimenti diventa impossibile andare avanti

Codice:

clear all
close all

tf=5; #tempo finale
ts=1000; #time steps
k=tf/ts; #passo temporale

#costanti
T0=5;
q0=0;
g=0; 
k0=100; %primo termine della polinomiale per la conduttività termica: k=k0+aT
k1=20;  %coefficiente della T

f=@(y) [y(2);-g./(k0 + k1*y(1))];
df=@(y) [1,0;0,(g*k1)/((k0+k1*y(1)).^2)];

#Eulero esplicito
y0=[0;2]; #DATO INIZIALE: corrisponde a u(0) e u'(0)=-g/k.

y=NaN(2,ts+1);
y(:,1)=y0;
for n=1:ts
   y(:,n+1)=y(:,n)+k*f(y(:,n));
end

t=linspace(0,tf,ts+1);
plot(t,y(1,:))


[luca.milano]

Ah quindi mi stai dicendo che anche la soluzione che ho trovato nel caso lineare fosse sbagliata? Eh si che sembrava tanto carina

[feddy]

Non so cosa debba venire, ma se non hai condizioni a un bordo non te le puoi certo inventare.

[luca.milano]

Hai ragione, probabilmente ho bisogno di capire prima un po' meglio la teoria e cercherò di leggere bene il testo che mi hai passato. A proposito, se hai altri testi, in particolare con esempi di codici matlab, sono ben accetti.

Non so cosa debba venire, ma se non hai condizioni a un bordo non te le puoi certo inventare

Non capisco proprio a cosa ti riferisci quando dici che io mi sia inventato una condizione a un bordo: io ho impostato la funzione al nodo $1$ e la funzione al nodo $2$ (nota la derivata prima in $0$) e la soluzione trovata nel caso lineare viene corretta (l'ho confrontata con la soluzione analitica e con altri risolutori, quindi sono sicuro che sia corretta).

[feddy]

Un testo recente è "Modellistica numerica per problemi differenziali - Quarteroni".
Per quanto riguarda l'ultima cosa: quello che intendo è che hai risolto un problema ai valori iniziali con un metodo per problemi ai limiti, visto che non conoscevi la condizione al bordo destro. Questo è sbagliato, che venga "giusta" poi è un altro conto. Oltretutto, quello che hai fatto te è stato mettere la condizione di Dirichlet in $0$ e nel nodo successivo quella di Neumann... ma converrai che così non hai risolto veramente il problema assegnato, visto che non hai modificato lo stesso nodo. Insomma, a scanso di equivoci: per ODE si usano metodi per ODE, per BVP differenze finite