欧拉近似方法求常微分方程.doc

欧拉近似方法求常微分方程
朱翼

1、编程实现以下科学计算算法，并举一例应用之。
“欧拉近似方法求常微分方程”
算法说明：
欧拉法是简单有效的常微分方程数值解法，欧拉法有多种形式的算法，其中简单欧拉法是一种单步递推算法。其基本原理为
对简单的一阶方程的初值问题：
y’=f(x,y)
其中 y(x0 )=y0
欧拉法等同于将函数微分转换为数值微分，由欧拉公式可得
yn+1 =y n+hf(x n ,y n)

程序代码：
function [tout,yout]=myeuler(ypfun,t0,tfinal,y0,tol,trace)
%初始化
pow=1/3;
if nargin<5,tol=1.e-3;end
if nargin<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.';
if trace %绘图
clc,t,h,y
end
while (th>t) %主循环
if t+h>tfinal,h=tfinal-t;end
% Compute the slopes
f=feval(ypfun,t,y);f=f(:);
%估计误差并设定可接受误差
delta=norm(h*f,'inf');
tau=tol*max(norm(y,'inf'),1.0);
%当误差可接受时重写解
if delta<=tau
t=t+h;
y=y+h*f;
k=k+1;
if k>length(tout)
tout=[tout;zeros(chunk,1)];
yout=[yout;zeros(chunk,length(y))];
end
tout(k)=t;
yout(k,:)=y.';
end
if trace
home,t,h,y
end
% Update the step size
if delta~=0.0
h=min(hmax,0.9*h*(tau/delta)^pow);
end
end
if (tunk=128;
tout=zeros(chunk,1);yout=zeros(chunk,length(y)); k=1;tout(k)=t;yout(k,:)=y.';

N

Y

trace= 1

N
输出t,h,y

N

tt

Y

Y

t+h>tfinal

h=tfinal-t

N

f=feval(ypfun,t,y); f=f(:);
delta=norm(h*f,'inf'); tau=tol*max(norm(y,'inf')1.0);

tau=tol*max(norm(y,'inf')1.0);

Y

delta <=tau

N
t=t+h; y=y+h*f; k=k+l;

k>length(tout)

tout=[tout;zeros(chunk,1)]; yout=[yout;zeros(chunk,length(y))];

tout(k)=t; yout(k,:)=y,';

Y

trace=1

输出home ,t ,h,y

N

Y

delta~=0.0

h=min(hmax,0.9*h*(tau/delta)^pow))

N

Y
t到其精确解为y（x）=x+exp(-x)。最后在同一图窗中分别画出两图。
程序代码：
f．m
function f=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('欧拉法','精确解')

例题流程图：

f=f(x,y) f=-y+x+1;

开始

调用函数myeuler ， [x,y]=myeuler(‘f’,0,1,1);
求出结果

y1=x+exp(-x)

调用函数plot 绘图，比较

结束

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 *

Columns 1 through 9

-0.0025 -0.0050 -0.0074 -0.0098 -0.0122 -0.0146 -0.0170 -0.0193 -0.0216

Columns 10 through 18

-0.0239 -0.0261 -0.0283 -0.0305 -0.0327 -0.0349 -0.0370 -0.0391 -0.0412

Columns 19 through 27

-0.0432 -0.0452 -0.0472 -0.0492 -0.0512 -0.0531 -0.0550 -0.0569 -0.0587

Columns 28 through 36

-0.0605 -0.0623 -0.0641 -0.0659 -0.0676 -0.0693 -0.0710 -0.0726 -0.0742

Columns 37 through 45

-0.0759 -0.0774 -0.0790 -0.0805 -0.0820 -0.0835 -0.0849 -0.0863 -0.0877

Columns 46 through 54

-0.0891 -0.0905 -0.0918 -0.0931 -0.0944 -0.0956 -0.0968 -0.0980 -0.0992

Columns 55 through 63

-0.1004 -0.1015 -0.1026 -0.1037 -0.1047 -0.1057 -0.1067 -0.1077 -0.1087

Columns 64 through 72

-0.1096 -0.1105 -0.1114 -0.1122 -0.1130 -0.1138 -0.1146 -0.1154 -0.1161

Columns 73 through 81

-0.1168 -0.1175 -0.1181 -0.1187 -0.1194 -0.1199 -0.1205 -0.1210 -0.1215

Columns 82 through 90

-0.1220 -0.1224 -0.1229 -0.1232 -0.1236 -0.1240 -0.1243 -0.1246 -0.1249

Columns 91 through 99

-0.1251 -0.1253 -0.1255 -0.1257 -0.1259 -0.1260 -0.1261 -0.1262 -0.1262

Columns 100 through 101

-0.1262 -0.1262

Y =

1.0e-004 *

Columns 1 through 9

-0.0002 -0.0007 -0.0015 -0.0025 -0.0037 -0.0052 -00069 -0.0088 -0.0109

Columns 10 through 18

-0.0133 -0.0159 -0.0188 -0.0218 -0.0251 -0.0286 -0.0323 -0.0362 -0.0403

Columns 19 through 27

-0.0446 -0.0492 -0.0539 -0.0588 -0.0639 -0.0692 -0.0747 -0.0804 -0.0863

Columns 28 through 36

-0.0924 -0.0986 -0.1050 -0.1116 -0.1184 -0.1253 -0.1324 -0.1397 -0.1471

Columns 37 through 45

-0.1547 -0.1624 -0.1703 -0.1783 -0.1865 -0.1949 -0.2034 -0.2120 -0.2208

Columns 46 through 54

-0.2297 -0.2388 -0.2479 -0.2572 -0.2667 -02762 -0.2859 -0.2957 -0.3057

Columns 55 through 63

-0.3157 -0.3258 -0.3361 -0.3465 -0.3569 -0.3675 -0.3782 -0.3890 -0.3998

Columns 64 through 72

-0.4108 -0.4218 -0.4330 -0.4442 -0.4555 -0.4669 -0.4784 -0.4899 -0.5015

Columns 73 through 81

-0.5132 -0.5249 -0.5367 -0.5486 -0.5606 -0.5726 -0.5846 -0.5967 -0.6088

Columns 82 through 90

-0.6210 -0.6333 -0.6456 -0.6579 -0.6703 -0.6827 -0.6951 -0.7075 -0.7200

Columns 91 through 99

-0.7325 -0.7451 -0.7576 -0.7702 -07828 -0.7954 -0.8080 -0.8206 -0.8332

Columns 100 through 101

-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 % 绘弯矩图

流程图

开始

L=1; P=1000; L1=1;
E=200*10^9;I=2*10^-5

调用linspace函数，将L分成100段

n1=L1/dx+1; M1 = -P*( L1-x(1:n1)); M2 = zeros(1,101-n1); M= [M1,M2];

调用cumsum函数求
A=cumsum(M)*dx/(E*I) Y=cumsum(A)*dx

结束

（2）计算积分
①
解题思路：exp（-x^2）是不可积函数，matlab中int积分显示不出积分结果，用到vpa函数控制其结果精度表示出来。
程序：
>> syms x
>> 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)。
程序：
>> syms a r
>> 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