這是非穩態一維熱傳導的方法,也叫古典顯格式。
如果是做數學建模,就別用了,這種方法計算量比較大,算的很慢,而且收斂不好。
但是如果實在沒辦法也能湊合用。
該改的地方我都用???代替了。
給個詳細解釋https://wenku.baidu.com/view/78a359d43b3567ec112d8a77.html?qq-pf-to=pcqq.group
1 function rechuandao() % 2 3 Llist = ??? 4 N = 1000; % 空間點數 5 M = 10000; 6 alfa = ??? % 導熱 / (密度*比熱) 7 lambda = 0.4; % 穩定條件 8 %************ 9 zone=Llist(1)/N;%空間步長 10 z=0:zone:Llist(1);%空間點坐標數組 11 z=z';%矩陣轉置 12 13 n = lambda * (Llist(1)/N)^2 / alfa(1); 14 15 TM = M * n; % 總時間 16 17 t=0:n:TM; %時間點坐標數組 18 t=t'; %矩陣轉置 19 20 %計算初值和邊值 21 T=zeros(N+1,M+1); %賦溫度矩陣初值 22 Ti=init_fun(z); %各空間點處的初始條件 23 To=border_funo(t); %內邊界條件 24 Te=border_fune(t); %外邊界條件 25 T(:,1)=Ti; %(初始條件 溫度都為37) 26 T(1,:)=To; %(邊界條件x=0處溫度為) 27 T(N+1,:)=Te; %(邊界條件x=L 處溫度為) 28 29 %用差分法求出溫度T與桿長L、時間t的關系 30 for k=1:M %時間點數 31 n=2; 32 while n<=N %空間點數 33 T(n,k+1)=lambda*(T(n+1,k)+T(n-1,k))+(-2*lambda+1)*T(n,k); 34 n=n+1; 35 end 36 end 37 %設置立體網格 38 for i=1:M+1 39 X(:,i)=z; 40 end 41 for j=1:N+1 42 Y(j,:)=t; 43 end 44 mesh(X,Y,T); 45 view([1 -1 1]); 46 xlabel('長度'); 47 ylabel('時間'); 48 zlabel('溫度'); 49 csvwrite('csv.csv',T) 50 51 function y=init_fun(z)%初值條件 52 y=??? 53 return 54 55 function y=border_funo(t)%z=0的邊界條件 56 y=??? 57 return 58 59 function y=border_fune(t)%z=L的邊界條件 60 y= ??? 61 return 62 %