声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2836|回复: 2

请教MATLAB语言编写关于欧拉方程的叠带程序,

[复制链接]
发表于 2005-11-24 23:08 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
   请教MATLAB语言编写关于欧拉方程的叠带程序,<BR>Un+1=<U>Un+hf(Tn,Un)</U>
回复
分享到:

使用道具 举报

发表于 2005-11-25 00:32 | 显示全部楼层

回复:(fenger)请教MATLAB语言编写关于欧拉方程的叠...

<P><br>function diff_numeric_solve_euler_1<br>%method1:euler法---------y(n+1)=yn+h*f(xn,yn)---&gt;显式计算公式<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))---&gt;隐式计算公式<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编辑过]

发表于 2009-1-16 11:58 | 显示全部楼层

请教MATLAB语言编写关于隐式欧拉方程的迭代程序

请教MATLAB语言编写关于隐式欧拉方程的迭代程序,而且这里的,,Ui和Bi都不是向量,而是矩阵的,Ui=U(i-1)+△t(A*Ui*Bi-Ui*Bi*inv(Ui)*A*Ui)   和Bi=Ui-B(i-1)*inv(U(i-1))  ,△t=0.01, A是已知的
急求高人!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-9-22 21:13 , Processed in 0.066049 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表