欧拉近似方法求常微分方程.docx
- 文档编号:2866796
- 上传时间:2023-05-04
- 格式:DOCX
- 页数:20
- 大小:238.96KB
欧拉近似方法求常微分方程.docx
《欧拉近似方法求常微分方程.docx》由会员分享,可在线阅读,更多相关《欧拉近似方法求常微分方程.docx(20页珍藏版)》请在冰点文库上搜索。
欧拉近似方法求常微分方程
欧拉近似方法求常微分方程
朱翼
1、编程实现以下科学计算算法,并举一例应用之。
“欧拉近似方法求常微分方程”
算法说明:
欧拉法是简单有效的常微分方程数值解法,欧拉法有多种形式的算法,其中简单欧拉法是一种单步递推算法。
其基本原理为
对简单的一阶方程的初值问题:
y’=f(x,y)
其中y(x0)=y0
欧拉法等同于将函数微分转换为数值微分,由欧拉公式可得
yn+1=yn+hf(xn,yn)
程序代码:
function[tout,yout]=myeuler(ypfun,t0,tfinal,y0,tol,trace)
%初始化
pow=1/3;
ifnargin<5,tol=1.e-3;end
ifnargin<6,trace=0;end
t=t0;
hmax=(tfinal-t)/16;
h=hmax/8;
y=y0(:
);
chunk=128;
tout=zeros(chunk,1);
yout=zeros(chunk,length(y));
k=1;
tout(k)=t;
yout(k,:
)=y.';
iftrace%绘图
clc,t,h,y
end
while(t
ift+h>tfinal,h=tfinal-t;end
%Computetheslopes
f=feval(ypfun,t,y);f=f(:
);
%估计误差并设定可接受误差
delta=norm(h*f,'inf');
tau=tol*max(norm(y,'inf'),1.0);
%当误差可接受时重写解
ifdelta<=tau
t=t+h;
y=y+h*f;
k=k+1;
ifk>length(tout)
tout=[tout;zeros(chunk,1)];
yout=[yout;zeros(chunk,length(y))];
end
tout(k)=t;
yout(k,:
)=y.';
end
iftrace
home,t,h,y
end
%Updatethestepsize
ifdelta~=0.0
h=min(hmax,0.9*h*(tau/delta)^pow);
end
end
if(t dish('Singularitylikely.') t end tout=tout(1: k); yout=yout(1: k,: ); 流程图: N Y N Y N N Y Y N Y N Y N Y N Y N 举例: 用欧拉法求y’=-y+x+1,y(0)=1。 解题思路: 首先建立例子所给函数的函数文件,调用函数myeuler,利用程序求解方程。 将欧拉法解和精确解比较,由方程f=-y+x+1可得到其精确解为y(x)=x+exp(-x)。 最后在同一图窗中分别画出两图。 程序代码: f.m functionf=f(x,y) f=-y+x+1; >>[x,y]=myeuler('f',0,1,1);%利用程序求解方程 >>y1=x+exp(-x);%方程f=-y+x+1的精确解 >>plot(x,y,'-b',x,y1,'-r')%在同一图窗将欧拉法解和精确解的图画出 >>legend('欧拉法','精确解') 例题流程图: 2、编程解决以下科学计算问题 (1) 解题思路: ◆建模: 材料力学中从弯矩求转角要经过一次不定积分,而从转角求挠度又要经过一次不定积分,在MATLAB中这却是非常简单的问题.可用cumsum函数作近似的不定积分,还可用更精确的函数cumtrapz来做不定积分。 本题用cumsum函数来做.解题的关键还是在于正确地列写弯矩方程。 本例中弯矩为 程序代码: >>L=1;P=1000;L1=1;%给出常数 E=200*10^9;I=2*10^-5; x=linspace(0,L,101);dx=L/100;%将x分100段 n1=L1/dx+1;%确定x=L1处对应的下标 M1=-P*(L1-x(1: n1));%第一段弯矩赋值 M2=zeros(1,101-n1);%第二段弯矩赋值(全为零) M=[M1,M2];%全梁的弯矩 A=cumsum(M)*dx/(E*I)%对弯矩积分求转角 Y=cumsum(A)*dx%对转角积分求挠度 A= 1.0e-003* Columns1through9 -0.0025-0.0050-0.0074-0.0098-0.0122-0.0146-0.0170-0.0193-0.0216 Columns10through18 -0.0239-0.0261-0.0283-0.0305-0.0327-0.0349-0.0370-0.0391-0.0412 Columns19through27 -0.0432-0.0452-0.0472-0.0492-0.0512-0.0531-0.0550-0.0569-0.0587 Columns28through36 -0.0605-0.0623-0.0641-0.0659-0.0676-0.0693-0.0710-0.0726-0.0742 Columns37through45 -0.0759-0.0774-0.0790-0.0805-0.0820-0.0835-0.0849-0.0863-0.0877 Columns46through54 -0.0891-0.0905-0.0918-0.0931-0.0944-0.0956-0.0968-0.0980-0.0992 Columns55through63 -0.1004-0.1015-0.1026-0.1037-0.1047-0.1057-0.1067-0.1077-0.1087 Columns64through72 -0.1096-0.1105-0.1114-0.1122-0.1130-0.1138-0.1146-0.1154-0.1161 Columns73through81 -0.1168-0.1175-0.1181-0.1187-0.1194-0.1199-0.1205-0.1210-0.1215 Columns82through90 -0.1220-0.1224-0.1229-0.1232-0.1236-0.1240-0.1243-0.1246-0.1249 Columns91through99 -0.1251-0.1253-0.1255-0.1257-0.1259-0.1260-0.1261-0.1262-0.1262 Columns100through101 -0.1262-0.1262 Y= 1.0e-004* Columns1through9 -0.0002-0.0007-0.0015-0.0025-0.0037-0.0052-0.0069-0.0088-0.0109 Columns10through18 -0.0133-0.0159-0.0188-0.0218-0.0251-0.0286-0.0323-0.0362-0.0403 Columns19through27 -0.0446-0.0492-0.0539-0.0588-0.0639-0.0692-0.0747-0.0804-0.0863 Columns28through36 -0.0924-0.0986-0.1050-0.1116-0.1184-0.1253-0.1324-0.1397-0.1471 Columns37through45 -0.1547-0.1624-0.1703-0.1783-0.1865-0.1949-0.2034-0.2120-0.2208 Columns46through54 -0.2297-0.2388-0.2479-0.2572-0.2667-0.2762-0.2859-0.2957-0.3057 Columns55through63 -0.3157-0.3258-0.3361-0.3465-0.3569-0.3675-0.3782-0.3890-0.3998 Columns64through72 -0.4108-0.4218-0.4330-0.4442-0.4555-0.4669-0.4784-0.4899-0.5015 Columns73through81 -0.5132-0.5249-0.5367-0.5486-0.5606-0.5726-0.5846-0.5967-0.6088 Columns82through90 -0.6210-0.6333-0.6456-0.6579-0.6703-0.6827-0.6951-0.7075-0.7200 Columns91through99 -0.7325-0.7451-0.7576-0.7702-0.7828-0.7954-0.8080-0.8206-0.8332 Columns100through101 -0.8459-0.8585 >>subplot(3,1,1),plot(x,M),grid%绘弯矩图 subplot(3,1,2),plot(x,A),grid%绘弯矩图 subplot(3,1,3),plot(x,Y),grid%绘弯矩图 流程图 结束 (2)计算积分 ① 解题思路: exp(-x^2)是不可积函数,matlab中int积分显示不出积分结果,用到vpa函数控制其结果精度,表示出来。 程序: >>symsx >>t=vpa(int(exp(-x.^2)./(1+x.^2),-inf,+inf),5) 结果: t=1.3433 ② 解题思路: 先建立内置函数,然后调用quad函数求积分。 程序: >>fun=@(x)tan(x)./x.^(0.7); >>quad('fun',0,1) 结果: ans=0.9063 ③ 解题思路: 先建立内置函数,然后调用quad函数求积分。 程序: >>fun=inline('exp(x)./(1-x.^2).^0.5'); >>quad('fun',0,1) 结果: ans=3.1044 ④ 解题思路: 这是一个二重积分,一般的方法是把二重积分化为二次积分,但由于二次积分命令 int(int(1+x+y.^2,y,-sqrt(2*x-x.^2),sqrt(2*x-x.^2)),x,0,1) 无法求出结果,故将其转化成极坐标。 积分函数为 r.*(1+r.*cos(a)+r.^2.*(sin(a)).^2) 其中a的范围: -pi/2到pi/2;r的范围: 0到2.*cos(a)。 程序: >>symsar >>int(int(r+r.^2.*cos(a)+r.^3.*(sin(a)).^2,r,0,2.*cos(a)),a,-pi/2,pi/2) 结果: ans=(9*pi)/4 (注: 本资料素材和资料部分来自网络,仅供参考。 请预览后才下载,期待您的好评与关注! )
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 近似 方法 微分方程