下面是matlab的程式碼，可以參考，其實主要部分是離散格式，掌握好網格劃分與差分格式，其他的都是程式設計問題了。%不知道你想比較什麼，直接執行出圖吧 N = 10;%建議取值10與30進行對比，取10時網格比較粗糙但效果可以，取30時結果發散。%希望不要追問隱式格式的，那個又是另一種了。T = 1; X = pi; dt = T/(N-1); dx = X/(N-1); t = 0:dt:T; x = 0:dx:X; co = dt/dx/dx; u = zeros(N); u(1,:) = sin(x); u(:,1) = 0; u(:,N) = 0; for i = 2:N for j = 2:N-1 u(i,j) = co*(u(i-1,j-1)+u(i-1,j+1)-2*u(i-1,j))+u(i-1,j); end end v = zeros(N); for i = 1:N for j = 1:N v(i,j) = exp(-t(i))*sin(x(j)); end end figure(); surf(x,t,u); figure(); surf(x,t,v);
顯示格式的效果比較差，而且對時間步長有限制，不然解就比較差了。
下面是matlab的程式碼，可以參考，其實主要部分是離散格式，掌握好網格劃分與差分格式，其他的都是程式設計問題了。%不知道你想比較什麼，直接執行出圖吧 N = 10;%建議取值10與30進行對比，取10時網格比較粗糙但效果可以，取30時結果發散。%希望不要追問隱式格式的，那個又是另一種了。T = 1; X = pi; dt = T/(N-1); dx = X/(N-1); t = 0:dt:T; x = 0:dx:X; co = dt/dx/dx; u = zeros(N); u(1,:) = sin(x); u(:,1) = 0; u(:,N) = 0; for i = 2:N for j = 2:N-1 u(i,j) = co*(u(i-1,j-1)+u(i-1,j+1)-2*u(i-1,j))+u(i-1,j); end end v = zeros(N); for i = 1:N for j = 1:N v(i,j) = exp(-t(i))*sin(x(j)); end end figure(); surf(x,t,u); figure(); surf(x,t,v);