《偏微分方程数值解》 实验报告
实验一
常微分方程的数值解法
一、 实验内容 1. 分别用欧拉法和改进的欧拉法（预报—校正格式）和四阶龙格库塔法求 解初值问题。 （h=0.1）
y ′ = −y + x + 1 ,0 ≤ x ≤ 1 (1) y 0 =1
并与真解作比较（列出表格，并画图） 2. 选取一种理论上收敛但是不稳定的算法对问题（1）进行计算，并与真解 作比较（列表，画图） 。 二、 实验结果分析
二 实验结果分析
通过这次实验我掌握了将得到的解进一步精确， 而且要学会比较这几种方法 的精确性，显然，四阶龙格库达比改进欧拉发精确，改进欧拉发比欧拉法精确。
程序结果：
0 1 0.0250000000000000 1.00030990808236 0.0500000000000000 1.00122941879688 0.0750000000000000 1.00274348398340 0.100000000000000 0.125000000000000 0.150000000000000 0.175000000000000 0.200000000000000 0.225000000000000 0.250000000000000 0.275000000000000 0.300000000000000 0.325000000000000 0.350000000000000 0.375000000000000 0.400000000000000 0.425000000000000 0.450000000000000 0.475000000000000 0.500000000000000 0.525000000000000 1.00483741833333 1.00749689930416 1.01070797154688 1.01445701892311 1.01873075361613 1.02351621605355 1.02880077891340 1.03457212180418 1.04081822141212 1.04752735143119 1.05468808618799 1.06228927773127 1.07032004691684 1.07876978334422 1.08762814863660 1.09688505571054 1.10653066070931 1.11655536294553
Euler('fun',0,1,1,10) 其中fun.m为
function z=fun(x,y) z=-y+x+1;
程序结果：
2、改进的欧拉法
编写 TranEuler.m 函数，代码如下
function TranEuler(fun,x0,y0,xn,N) x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; h=(xn-x0)/N; 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'] y1 = dsolve('Dy=-y+x+1','y(0)=1','x') x1 = 0:0.1:1; y2=x1 + exp(-x1); A=[x1',y2'] plot(x',y','O',x1',y2');
0.550000000000000 0.575000000000000 0.600000000000000 0.625000000000000 0.650000000000000 0.675000000000000 0.700000000000000 0.725000000000000 0.750000000000000 0.775000000000000 0.800000000000000 0.825000000000000 0.850000000000000 0.875000000000000 0.900000000000000 0.925000000000000 0.950000000000000 0.975000000000000 1
xlabel('x'); ylabel('y(x)'); title('Õæʵ½âÓë¸Ä½øµÄÅ·À-·¨Çó½â¶Ô±È'); grid on; legend('¸Ä½øÅ·À-½â','Õæʵ½â');程Leabharlann 结果：3、4 阶 Rk 法
编写 4 阶 Rk 函数，代码如下
clear all; clc; tic fun=inline('-y+x+1','x','y'); [x,y]=ode45(fun,[0,1],1); B=[x,y] y1 = dsolve('Dy=-y+x+1','y(0)=1','x') x1 = 0:0.1:1; y2=x1 + exp(-x1); A=[x1',y2'] plot(x,y,'O',x1',y2'); xlabel('x'); ylabel('y(x)'); title('Õæʵ½âÓë4½×RK·¨Çó½â¶Ô±È'); grid on; legend('4½×RK·¨½â','Õæʵ½â'); toc
1、欧拉法
第一步：编写 Euler.m 函数，代码如下
function Euler(fun,x0,y0,xn,N) x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; h=(xn-x0)/N; 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'] y1 = dsolve('Dy=-y+x+1','y(0)=1','x') x1 = 0:0.1:1; y2=x1 + exp(-x1); A=[x1',y2'] plot(x',y','O',x1',y2'); xlabel('x'); ylabel('y(x)'); title('Õæʵ½âÓëÅ·À-·¨Çó½â¶Ô±È'); grid on; legend('Å·À-½â','Õæʵ½â');
1.12694980786900 1.13770486830921 1.14881163717622 1.16026142740887 1.17204577466010 1.18415642032450 1.19658530493382 1.20932456811006 1.22236655099527 1.23570378089426 1.24932896529859 1.26323499184411 1.27741493050958 1.29186201972077 1.30656966094317 1.32153141864356 1.33674102227941 1.35219235372536 1.36787944238047