euler.m

关于matlab解热传导问题的解法

function diff_numeric_solve_euler_1
%method1:euler法---------y(n+1)=yn+h*f(xn,yn)--->显式计算公式
xn=0:.1:1;
y0=1; %赋初值
y_n_plus1=[];
for i=1:length(xn)-1
y_n_plus1=[y_n_plus1;y0+diff(xn([1:2]))*(y0-2*xn(i)/y0)];
y0=y_n_plus1(end);
end
y_n_plus1;
%method2:后退euler法---------y(n+1)=yn+h*f(x(n+1),y(n+1))--->隐式计算公式
xn1=0:.1:1;
y01=1; %赋初值
y_n2=[];
for i=1:length(xn1)-1
y02=y01+diff(xn1([1:2]))*(y01-2*xn1(i)/y01);
y_n2=[y_n2;y01+diff(xn1([1:2]))*(y02-2*xn1(i+1)/y02)];
y01=y_n2(end);
end
%上述两种方法的局部截断误差分别约为h^2*y''(xn)/2,-h^2*y''(xn)/2
%利用中心差分公式的梯形公式可以得到更好的精度
xn2=0:.1:1;
y11=1; 
y_n3=[];
for i=1:length(xn2)-1
y03=y11+diff(xn2([1:2]))*(y11-2*xn2(i)/y11);
y_n3=[y_n3;y11+.5*diff(xn2([1:2]))*(y11-2*xn2(i)/y11+y03-2*xn2(i+1)/y03)];
y11=y_n3(end);
end
%建立预测-校正系统的改进euler公式
xnp=0:.1:1;
y21=1; 
y_n4=[];
for i=1:length(xnp)-1
y04=y21+diff(xnp([1:2]))*(y21-2*xnp(i)/y21);
y_n4=[y_n4;y21+.5*diff(xnp([1:2]))*(y21-2*xnp(i)/y21+y04-2*xnp(i+1)/y04)];
y21=y_n4(end);
end
%精确的解析解
dy=dsolve('Dy=y-2*x/y','y(0)=1','x');
y_exact=subs(dy,'x',{.1:.1:1});
%ode45的解
f=inline('x-2*t/x','t','x');
[t,Y]=ode45(f,[0,1],1);
cc=flipdim(Y([end:-4:1]),1)';
y_ode=cc([2:end]);
%精度比较
y_compare=[y_n_plus1';y_n2';y_n3';y_n4';y_ode;y_exact];
y_compare_minor=[y_compare(6,:)-y_compare(1,:);...
y_compare(6,:)-y_compare(2,:);y_compare(6,:)-y_compare(3,:);...
y_compare(6,:)-y_compare(4,:);y_compare(6,:)-y_compare(5,:)];
var_y=var(y_compare_minor');