x_message(['Integrazione problema n.3 con';... 'schema BTCS.';'Effetti di addensamento nodi ';... ' '; ... 'Ver 1.0 (28/1/2003)']) function [s]=strech(N,P,Q) AN=N-1; DETA=1./AN; TQI=1./tanh(Q); for L=1:N AL=L-1; ETA=AL*DETA; DUM=Q*(1-ETA); DUM=1.-tanh(DUM)*TQI; s(L)=P*ETA+(1.-P)*DUM; end, endfunction cambia=2; while cambia==2, txt=['lambda';'mu ']; sig=x_mdialog('Parametri equazione campione',txt,['1.0';'0.1']); lambda=evstr(sig(1)); mu=evstr(sig(2)); txt=['N. nodi';'dt ';'N. max iter';'Res. min']; sig=x_mdialog('Parametri integrazione ',txt,['11';'0.499e-2';'1000';'1e-5']); N=evstr(sig(1)); dt=evstr(sig(2)); itermax=evstr(sig(3)); Resmin=evstr(sig(4)); //N=101; dxi=1/(N-1); c=lambda*dt/dxi; s=mu*dt/dxi^2; tschema=['FTCS (esplicito, I ord. tempo, II ord. spazio)';... 'BTCS (implicito, I ord. tempo, II ord. spazio)';... 'CN (implicito, II ord. tempo, II ord. spazio)']; //nschema=x_choose(tschema,['Scegli lo schema alle differenze finite']); nschema=2; testo(1)=strcat( ['------------------------------------------------------------']); testo(2) =strcat( ['Riassunto dei parametri equazione ed integrazione:']); testo(3) =strcat( ['Parametri equazione -----> lambda:' string(msprintf('%2.5f',lambda ))']); testo(4) =strcat( [' mu:' string(msprintf('%2.5f',mu ))']); testo(5) =strcat( [' mu/lambda:' string(msprintf('%2.5f',mu/lambda))']); testo(6) =strcat( ['Parametri integrazione --> N. nodi:' string(msprintf('%2.5f',N ))']); testo(7) =strcat( [' dx:' string(msprintf('%2.5f',dxi ))']); testo(8) =strcat( [' dt:' string(msprintf('%2.5f',dt ))']); testo(9) =strcat( [' c:' string(msprintf('%2.5f',c ))']); testo(10) =strcat( [' s:' string(msprintf('%2.5f',s ))']); testo(11)=strcat( [' N. iter max:' string(msprintf('%2.5f',itermax))']); testo(12)=strcat( [' Res. min.:' string(msprintf('%2.5f',Resmin))']); testo(13)=strcat( [' Schema di integrazione:' string(msprintf('%s',tschema(nschema) ))']); testo(14)=strcat( ['------------------------------------------------------------']); testo(15)=strcat( ['Vuoi cambiare i parametri?']); cambia= x_message(testo,['No' 'Yes']); end, cambia=2; while cambia==2, txt=['P ([1.01 ..1.99] '] sig=x_mdialog('Parametro di addensamento ',txt,['1.01']); P=evstr(sig(1)); //P=1.8; Q=2.0; X=strech(N,P,Q); Xi=[0 : dxi : 1.]'; gn=X; gnp1=zeros(1,N)'; for i=2:N-1 dum= (X(i+1)-X(i-1))/(2*dxi); dxidx(i)=1./dum; end dxidx(1)=3*dxidx(2)-3*dxidx(3)+dxidx(4); dxidx(N)=3*dxidx(N-1)-3*dxidx(N-2)+dxidx(N-3); for i=2:N-1 // dxidx2(i)=(dxidx(i+1)-dxidx(i-1))/(X(i+1)-X(i-1)); dxidx2(i)=(dxidx(i+1)-dxidx(i-1))/(Xi(i+1)-Xi(i-1))*dxidx(i); end //dxidx=1.; //dxidx2=0.; iter=2000; //if nschema==1 then //FTCS // A=zeros(N,N); // A(1,1)=1.; // for i=2:N-1 // d=mu*dt/dxi*dxidx2(i); // A(i,i+1)=s*dxidx(i)^2-c*dxidx(i)/2+d/2.; // A(i,i)=1-2*s*dxidx(i)^2; // A(i,i-1)=(c*dxidx(i)/2+s*dxidx(i)^2-d/2.); // end // A(N,N)=1.; // L=A; //elseif nschema==2 then // //BTCS A=zeros(N,N); A(1,1)=1.; for i=2:N-1 d=mu*dt/dxi*dxidx2(i); A(i,i+1)=c*dxidx(i)/2-s*dxidx(i)^2-d/2.; A(i,i)=1+2*s*dxidx(i)^2; A(i,i-1)=-(c*dxidx(i)/2+s*dxidx(i)^2)+d/2; end A(N,N)=1.; L=inv(A); //else // //CN // A=zeros(N,N); // B=zeros(N,N); // A(1,1)=1.; // B(1,1)=1.; // for i=2:N-1 // A(i,i+1)=1/2*(c/2-s); // B(i,i+1)=1/2*(s-c/2); // A(i,i)=1+s; // B(i,i)=1-s; // A(i,i-1)=-1/2*(c/2+s); // B(i,i-1)=1/2*(c/2+s); // end // A(N,N)=1.; // B(N,N)=1.; // L=inv(A)*B; //end, xset("window",0) xtitle('Evoluzione della soluzione ', 'x', 'g'); //xset("window",1) //xtitle('Storia del residuo ', 'iter.', 'Res.'); xset("window",0) xselect(); kp=xget("pixmap");xset("pixmap",1); xset("pixmap",kp); kp=xget("pixmap");xset("pixmap",1); xselect();plot2d(X,gn,1,rect=[0.,0.,1.,1.]); xset("mark",-1,2); t=0.; res=1e+39; nn=1; gnp1max=0.; t=0.; while nn Resmin, gnp1=L*gn; gnp1max=max(gnp1); if gnp1max > 1.01 then x_message(['Attenzione g>1' ; 'Probabile schema di integrazione instabile ']); break; end; xset("window",0) xset("wwpc"); xtape("clear",0); xtitle('Evoluzione della soluzione ', 'x', 'g'); plot2d(X,gnp1,-1,rect=[0.,0.,1.,1.]); xset("wshow"); // xset("window",1); res=sum(abs((gnp1-gn))/dt)/N; // plot2d(nn,res,-1,logflag="nl", rect=[0,Resmin,itermax,10.]); t=t+dt; gn=gnp1; nn=nn+1; end, xset("window",0); ge=(exp(lambda*X/mu)-1.0)/(exp(lambda/mu)-1.0); plot2d(X,ge,2,rect=[0.,0.,1.,1.]); xset("wshow"); xset("window",1); //xtape("replay",1); //errore numerico xset("window",2); xtitle('Errore globlale', 'x', 'e'); plot2d(X,ge-gnp1); //errore di troncamento xset("window",1); xtitle('Errore di troncamento', 'x', 'Et'); et=(ge-L*ge)/dt; plot2d(X,et); //cont1=x_dialog('Lettura soluzione in un punto?','Y'); //if cont1=="Y" then // cont2="Y"; // while cont2=="Y", // inode=evstr(x_dialog('Scegli il nodo da valutare(1-N)','0')); // testo1(1) =strcat( ['Soluzione numerica e soluzione esatta a']); // testo1(2) =strcat( [' t: ' string(msprintf('%2.5f', t ))']); // testo1(3) =strcat( [' X: ' string(msprintf('%2.5f', X(inode) ))']); // testo1(4) =strcat( ['----------------------------------------------------------------']); // testo1(5) =strcat( ['soluzione --> Esatta:' string(msprintf('%2.5f',ge(inode) ))']); // testo1(6) =strcat( ['----------------------------------------------------------------']); // testo1(7) =strcat( ['Soluzione -->Numerica:' string(msprintf('%2.5f',gnp1(inode) ))']); // testo1(8) =strcat( ['Ottenuta con:' string(msprintf('%s',tschema(nschema))) ]); // testo1(9) =strcat( [' dx:' string(msprintf('%f',dx )) ' dt:' string(msprintf('%f',dt ))]); // x_message(testo1); // cont2=x_dialog('Altra valutazione ?','Y'); // end, //end, //xdel(0); //xdel(1); //xdel(2); //quit; testo1(1)='cambiare addensamento'; cambia= x_message(testo1,['No' 'Yes']); end, xdel(0); //xdel(1); xdel(2); quit; //Res=1; //xselect(); //kp=xget("pixmap");xset("pixmap",1) //xset("pixmap",kp); //kp=xget("pixmap");xset("pixmap",1); //xselect();plot2d(X',U',0,rect=[0.,0.,1.,1.]); //while Res > 1.0e-6 // for i=2:N-1 // Un(i)=U(i)-CFL*(U(i)-U(i-1))+CFL/Rex*(U(i+1)-2*U(i)+U(i-1)); // end // Res=sum(abs(U-Un))/N/CFL*dx; // Res; // U=Un; // xset("wwpc"); // plot2d(X',U',0,rect=[0.,0.,1.,1.]) // xset("wshow"); //end //xset("pixmap",kp);