微分方程数值解课程设计题目
1（30分）分别用Euler法、改进的Euler法及四阶的龙格库塔法求解初值问题：
h分别取0.1和0.2，要求计算并绘出图形，最后比较三种算法的精度。
（1）.先构建初值问题的函数
M文件：
function z=fun(x,y)
z=y-2*x/y;
（2）.Euler法
M文件：
function E=Euler(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));
end
T=[x;y]
1 .当h=0.1时
>> Euler('fun',0,1,0.1,10)
T = Columns 1 through 10
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000
1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 1.6498 1.7178
Column 11
1.0000
1.7848
②.当h=0.2时
>> Euler('fun',0,1,0.2,5)
T = 0 0.2000 0.4000 0.6000 0.8000 1.0000
1.0000 1.2000 1.3733 1.5315 1.6811 1.8269
（3）.改进的Euler法：
M文件：Euler_modify.m
function E=Euler_modify(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
z0=y(n)+h*feval(fun,x(n),y(n));
y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));
end
T=[x;y]
1 .当h=0.1时
>> Euler_modify('fun',0,1,0.1,10)
T = Columns 1 through 9
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000
1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165
Columns 10 through 11
0.9000 1.0000
1.6782 1.7379
②.当h=0.2时
>> Euler_modify('fun',0,1,0.2,5)
T = 0 0.2000 0.4000 0.6000 0.8000 1.0000
1.0000 1.1867 1.3483 1.4937 1.6279 1.7542
（4）.四阶R_K方法：
M文件：
function [x,y]=Rk_N4(f,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(f,x(n),y(n));
k2=h*feval(f,x(n)+1/2*h,y(n)+1/2*k1);
k3=h*feval(f,x(n)+1/2*h,y(n)+1/2*k2);
k4=h*feval(f,x(n)+h,y(n)+k3);
y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
end
1 .当h=0.1时
>> [x,y]=Rk_N4('fun',0,1,0.1,10)
x = Columns 1 through 9
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000
Columns 10 through 11
0.9000 1.0000
y = Columns 1 through 9
1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125
Columns 10 through 11
1.6733 1.7321
②.当h=0.2时
>> [x,y]=Rk_N4('fun',0,1,0.2,5)
x = 0 0.2000 0.4000 0.6000 0.8000 1.0000
y = 1.0000 1.1832 1.3417 1.4833 1.6125 1.7321
（5）作图：
1): h=0.1
>> x=0:0.1:1;
>> y_euler=[1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 1.6498 1.7178 1.7848];
>> y_euler_modify=[1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379];
>> y_RK_N4=[1.0000 1.0954 1.1832 &n伏羲画八卦
bsp; 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321];
plot(x,y_euler,'bo:',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');ylabel('y轴');text(0.3,1.3,'Euler法');text(0.4,1.35,'Euler改进法');text(0.5,1.4,'R-k法'); text(0.8,1.02,'作者：李靖');text(0.8,1.07,'日期：2010.6.19');text(0.4,1.75,'h=0.1');legend('Euler','Euler改进法','R_K法');grid on
2): h=0.2
>> x=0:0.2:1;
>> y_euler=[1.0000 1.2000 1.3733 1.5315 1.6811 1.8269];
>> y_euler_modify=[1.0000 1.1867 1.3483 1.4937 1.6279 1.7542];
>> y_RK_N4=[ 1.0000 1.1832 1.3417 1.4833 1.6125 1.7321];
plot(x,y_euler,'bo:',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');ylabel('y轴');text(0.3,1.3,'Euler法');text(0.4,1.35,'Euler改进法');text(0.5,1.4,'R-k法'); text(0.8,1.02,'作者：xx');text(0.8,1.07,'日期：2010.6.19'); text(0.4,1.75,'h=0.2');legend('Euler','Euler改进法','R_K法');grid on
2.（20分）编写一个程序用tylor级数法求解问题：
取tylor级数法的截断误差为O()，即要用u(t)，u’(t)，…，u(20)(t)的值。（提示：可用一个简单的递推公式来获得u(n)(t)，n=1,2，……）。
（1）.先利用递推公式计算u（t）的各阶导数
M文件：
differ.m
function D=differ(N)
y=zeros(1,N+1);
syms x y;
y(1)=x*y;
for n=1:N
y(n+1)=y(n)*x+diff(y(n),x);
end
y
输入（由于随着导数的阶数增高，项越来越多。所以只举例计算到5阶）：
>> differ(5)
y =
[x*y, x^2*y+y, (x^2*y+y)*x+2*x*y, 冬笋怎么做好吃又简单
((x^2*y+y)*x+2*x*y)*x+3*x^2*y+3*y, (((x^2*y+y)*x+2*x*y)*x+3*x^2*y+3*y)*x+(3*x^2*y+3*y)*x+(x^2*y+y)*x+8*x*y, ((((x^2*y+y)*x+2*x*y)*x+3*x^2*y+3*y)*x+(3*x^2*y+3*y)*x+(x^2*y+y)*x+8*x*y)*x+((3*x^2*y+3*y)*x+(x^2*y+y)*x+8*x*y)*x+((x^2*y+y)*x+2*x*y)*x+15*x^2*y+15*y]
