Questa è l'equazione:
$ { ( (partial^2 T)/(partial x^2) + g/k=0 ),( T(0)=T_0 ),( T(L)=T_L ):} $
dove il termine $ k $ è pari a:
$ k=k0+a*T $
Vi riporto lo script principale:
- Codice:
L=1; %lunghezza
N=50; %numero nodi
dx=L/N; %intervallo spaziale
T0=0; %temperatura in x=0
TL=0; %temperatura in x=L
g=50; %termine di generazione
k0=10; %primo termine della polinomiale per la conduttività termica: k=k0+aT
a=10; %coefficiente della T
save indiretto
%METODO DI NEWTON
T=zeros(N,1);
x1=ones(N,1);
n=0;
while norm(T-x1)>1e-5
n=n+1;
x1=T;
r=jac2(T)\sist2(T);
T=T-r;
end
T;
n;
La function per il vettore delle equazioni è la seguente:
- Codice:
function F=sist2(T)
%devo creare il vettore colonna del sistema: F=0
load indiretto
F(1)=T(1)-T0;
F(N)=T(L)-TL;
for j=2:N-1
k(j)=k0+(10^-1)*T(j);
F(j)=T(j+1)-2*T(j)+T(j-1)+g/k(j)*dx^2;
end
F=F';
Lo Jacobiano è invece questo:
- Codice:
function J=jac2(T)
load indiretto
J=zeros(N,N);
J(1,1)=1; J(N,N)=1;
for j=2:N-1
k(j)=k0+10^-1*T(j);
J(j,j-1)=1;
J(j,j)=-2-a*g/(k0+a*T(j))^2;
J(j,j+1)=1;
end
Vedete qualche errore? Ci sono metodi migliori per risolvere sistemi non lineari?
.
