这是非稳态一维热传导的方法,也叫古典显格式。
如果是做数学建模,就别用了,这种方法计算量比较大,算的很慢,而且收敛不好。
但是如果实在没办法也能凑合用。
该改的地方我都用???代替了。
给个详细解释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 %