|
回复:(fenger)请教MATLAB语言编写关于欧拉方程的叠...
<P><br>function diff_numeric_solve_euler_1<br>%method1:euler法---------y(n+1)=yn+h*f(xn,yn)--->显式计算公式<br>xn=0:.1:1;<br>y0=1; %赋初值<br>y_n_plus1=[];<br>for i=1:length(xn)-1<br>y_n_plus1=[y_n_plus1;y0+diff(xn([1:2]))*(y0-2*xn(i)/y0)];<br>y0=y_n_plus1(end);<br>end<br>y_n_plus1;<br>%method2:后退euler法---------y(n+1)=yn+h*f(x(n+1),y(n+1))--->隐式计算公式<br>xn1=0:.1:1;<br>y01=1; %赋初值<br>y_n2=[];<br>for i=1:length(xn1)-1<br>y02=y01+diff(xn1([1:2]))*(y01-2*xn1(i)/y01);<br>y_n2=[y_n2;y01+diff(xn1([1:2]))*(y02-2*xn1(i+1)/y02)];<br>y01=y_n2(end);<br>end<br>%上述两种方法的局部截断误差分别约为h^2*y''(xn)/2,-h^2*y''(xn)/2<br>%利用中心差分公式的梯形公式可以得到更好的精度<br>xn2=0:.1:1;<br>y11=1; <br>y_n3=[];<br>for i=1:length(xn2)-1<br>y03=y11+diff(xn2([1:2]))*(y11-2*xn2(i)/y11);<br>y_n3=[y_n3;y11+.5*diff(xn2([1:2]))*(y11-2*xn2(i)/y11+y03-2*xn2(i+1)/y03)];<br>y11=y_n3(end);<br>end<br>%建立预测-校正系统的改进euler公式<br>xnp=0:.1:1;<br>y21=1; <br>y_n4=[];<br>for i=1:length(xnp)-1<br>y04=y21+diff(xnp([1:2]))*(y21-2*xnp(i)/y21);<br>y_n4=[y_n4;y21+.5*diff(xnp([1:2]))*(y21-2*xnp(i)/y21+y04-2*xnp(i+1)/y04)];<br>y21=y_n4(end);<br>end<br>%精确的解析解<br>dy=dsolve('Dy=y-2*x/y','y(0)=1','x');<br>y_exact=subs(dy,'x',{.1:.1:1});<br>%ode45的解<br>f=inline('x-2*t/x','t','x');<br>[t,Y]=ode45(f,[0,1],1);<br>cc=flipdim(Y([end:-4:1]),1)';<br>y_ode=cc([2:end]);<br>%精度比较<br>y_compare=[y_n_plus1';y_n2';y_n3';y_n4';y_ode;y_exact];<br>y_compare_minor=[y_compare(6,:)-y_compare(1,:);...<br>y_compare(6,:)-y_compare(2,:);y_compare(6,:)-y_compare(3,:);...<br>y_compare(6,:)-y_compare(4,:);y_compare(6,:)-y_compare(5,:)];<br>var_y=var(y_compare_minor');<br><br>转simwe上的一个例子</P>
[此贴子已经被作者于2005-11-25 0:48:13编辑过]
|
|