數值計算求解動態熱傳導方程

和求解靜態熱傳導方程不同的是，動態代表溫度會隨時間而變化，這就要用到之前學的幾個初值函式，先來複習一下這幾個初值函式

function [lhs,rhs]=ost(theta,timestep,m,b,c,sol)

lhs=m-theta*timestep*b(1)
rhs=(m+(1-theta))*timestep*b(2))*sol+timestep*(theta*c(1)+(1-theta))*c(2)
end%之前b和c都是陣列元素引用，但這裡是矩陣，b(1)b(2)在這裡都是b，同時動態熱傳導方程裡面的c是恒為0
ab2function [lhs,rhs] = ab2(timestep,m,b,c,sol)
lhs=m;
rhs=(m+3/2*timestep*b(1))*sol(1)+timestep/2*(3*c(1)-b(2)*sol(2)-c(2));
endam3
function [lhs,rhs] = am3(timestep,m,b,c,sol)
lhs=m-5/12*timestep*b(1);
rhs=(m+2/3*timestep*b(2))*sol(1)+timestep/12*(5*c(1)+8*c(2)-b(3)*sol(2)-c(3));
endbdf2
function [lhs,rhs] = bdf2(timestep,m,b,c,sol)
lhs=3/2*m-timestep*b;
rhs=2*m*sol(1)-1/2*m*sol(2)+timestep*c;
end

function val = linquadref(xi,eta)

n1=1./4*(1-xi)*(1-eta)
n2=1./4*(1+xi)*(1-eta)
n3=1./4*(1+xi)*(1+eta)
n4=1./4*(1-xi)*(1+eta)
val=[n1 n2 n3 n4]'
end

function [elemat,elevec]=evaluate_instat(elenodes,gpx,gpw,elesol,eleosol,timint_m,timestep,theta,firststep)

function [elemat,elevec]=evaluate_instat(elenodes,gpx,gpw,elesol,eleosol,timint_m,timestep,theta,firststep)

elemat=zeros(4,4);
lamda=48;
p=7800;
c=452;
m=zeros(4,4);
b=zeros(4,4);
for i=1:4
for j=1:4
s1=0;
s2=0;
for k=1:size(gpw)
xi=gpx(k,1);
eta=gpx(k,2);
deriv=linquadderivref(xi,eta);
val=linquadref(xi,eta);
[j,detj,invj]=getjacobian(elenodes,xi,eta);
nv=val(i,:);
nt=val(j,:);
gnv=deriv(i,:)*invj  %gnv gradient of n
gnt=deriv(j,:)*invj  
s1=s1+p*c*nv*nt*detj*gpw(k);
s2=s2+lamda*dot(gnv,gnt)*detj*gpw(k);
endm(i,j)=s1;
b(i,j)=-s2
% [lhs(i,j),rhs(i)] = ost(0.66,1000,m(i,j),[b(i,j),b(i,j)],[0,0],elesol(i))
endc=0;
if timit_m==1
if timint_m==1
[lhs,rhs] = ost(theta,timestep,m,b,c,elesol);
elseif timint_m==2
[lhs,rhs] = ab2(timestep,m,b,c,[elesol,eleosol]);
elseif timint_m==3
[lhs,rhs] = am3(timestep,m,b,c,[elesol,eleosol]);
elseif timint_m==4
[lhs,rhs] = bdf2(timestep,m,b,c,[elesol,eleosol]);
else
disp('false')
endelemat=lhs;
elevec=rhs;
end

clear;

r=0.02;
b=0.3;
h=0.3;
ts=5000;
t=0:500:5000;
eleosol=300*ones(18,1);
timint_m=1;
timestep=500;
theta=0.5;
firststep=1;
t(:,1)=300*ones(18,1);
elesol=300*ones(18,1);
nodes=[0 0;b/3 0;2*b/3 0;b 0;
0 h/3;b/3 h/3;2*b/3 h/3;b h/3;
0 2*h/3;b/3 2*h/3;2*b/3 2*h/3;b-r*sin(pi/6) h-r*cos(pi/6);
b h-r;b-r*cos(pi/6) h-r*sin(pi/6);0 h;b/3 h;
b/2 h;b-r h];
elements=[1 2 6 5;2 3 7 6;3 4 8 7;5 6 10 9;6 7 11 10;11 7 12 14;
7 8 13 12;9 10 16 15;10 11 17 16;11 14 18 17];
dbc=[1 600;2 600;3 600;4 600;
12 300;13 300;14 300;18 300];
for i=1:length(t)-1
sysmat=zeros(18,18);
rhs=zeros(18,1);
for n=1:size(elements,1)
[elemat,elevec]=evaluate_instat(nodes(elements(n,:),:),gx2dref(2),gw2dref(2),elesol(elements(n,:)),eleosol,timint_m,timestep,theta,firststep);
[sysmat,rhs]=assemble(elemat,elevec,sysmat,rhs,elements(n,:));
end[sysmat,rhs]=assigndbc(sysmat,rhs,dbc);
t(:,i+1)=sysmat\rhs;
elesol=t(:,i+1);
endquadplot(nodes,elements,t(:,length(t)));
shading interp;
grid on;
colormap(hot);
colobar

...

elesol=t(:,i+1);%**一直到這裡是一樣的
if t(15,i+1)>=450|t(16,i+1)>=450|t(17,i+1)>=450|t(18,i+1)>=450
break
endend
zeit=t(i+1);
quadplot(nodes,elements,t(:,i+1));
shading interp;
grid on;
colormap(hot);
colorbar

Matlab系列教程 數值計算 求方差和標準差

首先，什麼是方差和標準差？方差，是在概率論和統計方差衡量隨機變數或一組資料時離散程度的度量，統計中的方差 樣本方差 是每個樣本值與全體樣本值的平均數之差的平方值的平均數。在許多實際問題中，研究方差即偏離程度有著重要意義。標準差，中文環境中又常稱均方差，是離均差平方的算術平均數的平方根。標準差是方差的...

Matlab系列教程之數值計算求方差和標準差

首先，什麼是方差和標準差？方差，是在概率論和統計方差衡量隨機變數或一組資料時離散程度的度量，統計中的方差 樣本方差 是每個樣本值與全體樣本值的平均數之差的平方值的平均數。在許多實際問題中，研究方差即偏離程度有著重要意義。標準差，中文環境中又常稱均方差，是離均差平方的算術平均數的平方根。標準差是方差的...

數值計算 線性方程組求解（0）

本專題將講述以多種方式求解線性方程組，也作為本人在 數值計算與優化 課程中學到知識的總結與具體 實現。主要用到資料為 數值計算方法 第3版 解線性代數組是科學研究與工程計算中經常遇到的問題。此專題討論以下n階線性方程組 a 11x1 a12x2 a 1nxn b1a 21x1 a22x2 a 2nx...