求下列微分方程初值问题的数值解

1. 求下列微分方程初值问题的数值解：

(1)⎧dy 23⎪=1+x +y 在[0,1]上的解； ⎨dx ⎪y (0)=0⎩

⎧dy ⎪=-y +x +1⎨dx ⎪⎩y (0)=2(2)1≤x ≤1.6

解：(1) 建立一个m 文件euler .m ，（具体文件见附件文件夹）

function [x, y]=euler(dyfun,xspan,y0,N)

if nargin

N=100;

end

h=(xspan(2)-xspan(1))/N;

x = xspan(1):h:xspan(2);

y(1)=y0;

for n=1:(length(x)-1)

k1=feval(dyfun,x(n),y(n));

y(n+1)=y(n)+h*k1;

k2=feval(dyfun,x(n+1),y(n+1));

y(n+1)=y(n)+h*(k1+k2)/2;

end

x=x';y=y';

disp(' x点 x点处的微分值' );

for i = 1:length(x)

fprintf(' %7f %7f\n',x(i),y(i));

end

再建立一个m 文件dyfun.m(见附件文件夹)

function z = fdyfun(x,y)

z = x^2 + y^3;

在MATLAB 命令窗口调用euler.m 文件，这里取N=10

>> [x, y]=euler(@dyfun,[0,1],0,10);

执行结果是如下：

x点 x点处的微分值

0.000000 0.000000

0.100000 0.000500

0.200000 0.003000

0.300000 0.009500

0.400000 0.022000

0.500000 0.042504

0.600000 0.073023

0.700000 0.115607

0.800000 0.172408

0.900000 0.245829

1.000000 0.338842

（2）. 只要改动dyfun.m 文件就可以了

function z = fdyfun(x,y)

z = -y + x + 1;

然后取N=20，在MATLAB 命令窗口调用euler.m 文件；

>>[x, y]=euler(@dyfun,[1,1.6],2,20);

执行结果如下：

x 点 x 点处的微分值

1.000000 2.000000

1.030000 2.000450

1.060000 2.001773

1.090000 2.003944

1.120000 2.006937

1.150000 2.010728

1.180000 2.015293

1.210000 2.020610

1.240000 2.026657

1.270000 2.033411

1.300000 2.040852

1.330000 2.048960

1.360000 2.057715

1.390000 2.067097

1.420000 2.077089

1.450000 2.087672

1.480000 2.098829

1.510000 2.110543

1.540000 2.122797

1.570000 2.135575

1.600000 2.148862

2．分别用二分法、不动点迭代法、牛顿法求下列方程的根；

(1)

(2)x 3-x -1=0, x ∈[1,2],ε=10-10（用二分法） x 5-4x -2=0的根，x 0=1.5 （用不动点迭代）

（3）. 用Newton 法求x 3-x -1=0, x ∈[1,2],ε=10-10

解：（1）建立一个m 文件erfen.m, （见附件文件夹）

%二分法

function x=erfen(fun,a,b,ep)

%用途：用二分法求非线性方程f(x)=0有根区间[a,b]中的一个根;

%===============输入=============================%

%格式：x=erfen(fun,a,b,ep)

% fun为函数表达式，

% a, b为区间左右端点，

% ep为精度，

%================================================%

%===============输出=============================%

% x返回近似根

%================================================%

if nargin

error('输入至少有被根函数和根所在区间的上界和下界');

elseif nargin == 3

ep = 1e-15;

end

x=(a+b)/2.0;

i=1;

k=0;

m=abs(subs(fun,x))>ep;

n = (b-a>ep);

while m|n

f1 = subs(fun,a);

f2 = subs(fun,x);

if f1*f2

b=x;

else

a=x;

end

x=(a+b)/2.0;

k=k+1;

X(i)=x;

M(i) = k;

i=i+1;

m=abs(subs(fun,x))>ep;

n = (b-a>ep);

end

disp('迭代次数 对应次数得到的 X 的值');

for i=1:length(X)

fprintf(' %4d %-12f\n',M(i),X(i));

end

disp(['一共迭代的次数是：',num2str(k)]);

disp('最终的近似解是：');

然后建立一个m 文件fun.m （见附件文件夹）

function y = fun(x)

syms x

y = x^3-x-1;

在MATLAB 命令窗口调用erfen.m 文件

>> x=erfen(fun,1,2,1e-10)

执行结果如下：

迭代次数 对应次数得到的 X 的值

1 1.250000

2 1.375000

3 1.312500

4 1.343750

5 1.328125

6 1.320313

7 1.324219

8 1.326172

9 1.325195

10 1.324707

11 1.324951

12 1.324829

13 1.324768

14 1.324738

15 1.324722

16 1.324715

17 1.324718

18 1.324717

19 1.324718

20 1.324718

21 1.324718

22 1.324718

23 1.324718

24 1.324718

25 1.324718

26 1.324718

27 1.324718

28 1.324718

29 1.324718

30 1.324718

31 1.324718

32 1.324718

33 1.324718

34 1.324718

一共迭代的次数是：34

最终的近似解是：

x =

1.3247

（2）把x 5-2或ϕ(x ) =来进行迭代。 x -4x -2=0变为ϕ(x ) =

45

建立m 文件maiter.m(见附件文件夹)

%不动点迭代

function x=maiter(phi,x0,ep,N)

% 用途：用简单迭代法求非线性方程f(x)=0有根区间[a,b]中的一个根 % 格式：x= maiter(phi,x0,ep,N)

% fun为phi(x)的函数句柄,

% x0为初值,

% ep为精度(默认1e-4),

% N为最大迭代次数(默认500),

% x返回近似根

if nargin

if nargin

k=0;

while k

x=feval(phi,x0);

if abs(x-x0)

break;

end

x0=x;k=k+1;

end

if k==N, warning('已达迭代次数上限'); end

disp(['k=',num2str(k)])

再建立一个m 文件fun1.m(见附件文件夹)

function y = fun1(x)

y = (x^5-2)/4;

%y = (4*x+2)^(1/5);

在MATLAB 命令窗口调用maiter.m 文件,N 取10 x 5-2这里用ϕ(x ) =迭代 4

>> x=maiter(@fun1,1.5,10)

k=0

x =

1.3984

这里用ϕ(x ) =

>> x=maiter(@fun1,1.5,10)

k=0

x =

1.5157

（3）建立一个m 文件Newton.m(见附件文件夹)

%Newton迭代法

function x=newton(fun,dfun,x0,ep,N)

%用途：用牛顿法求解非线性方程f(x)=0

%格式：x=newton(fun,dfun,x0,ep,N) fun和dfun 分别为表示f(x)的 %函数句柄, x0为迭代初值, ep为精度(默认1e-4), N为最大迭代 %次数(默认为500), x返回近似根

if nargin

if nargin

k=0;

while k

x=x0-subs(fun,x0)/subs(dfun,x0);

if abs(x-x0)

break;

end

x0=x; k=k+1;

end

if k==N, warning('已达迭代次数上限'); end

disp(['k=',num2str(k)])

再建立一个m 文件dfun.m

function y1 = dfun(x)

syms x

y1 = diff(fun,x);

在MATLAB 命令窗口调用得

>> x=newton(fun,dfun,1,1e-10,20) k=5

x =

1.3247

1. 求下列微分方程初值问题的数值解：

(1)⎧dy 23⎪=1+x +y 在[0,1]上的解； ⎨dx ⎪y (0)=0⎩

⎧dy ⎪=-y +x +1⎨dx ⎪⎩y (0)=2(2)1≤x ≤1.6

解：(1) 建立一个m 文件euler .m ，（具体文件见附件文件夹）

function [x, y]=euler(dyfun,xspan,y0,N)

if nargin

N=100;

end

h=(xspan(2)-xspan(1))/N;

x = xspan(1):h:xspan(2);

y(1)=y0;

for n=1:(length(x)-1)

k1=feval(dyfun,x(n),y(n));

y(n+1)=y(n)+h*k1;

k2=feval(dyfun,x(n+1),y(n+1));

y(n+1)=y(n)+h*(k1+k2)/2;

end

x=x';y=y';

disp(' x点 x点处的微分值' );

for i = 1:length(x)

fprintf(' %7f %7f\n',x(i),y(i));

end

再建立一个m 文件dyfun.m(见附件文件夹)

function z = fdyfun(x,y)

z = x^2 + y^3;

在MATLAB 命令窗口调用euler.m 文件，这里取N=10

>> [x, y]=euler(@dyfun,[0,1],0,10);

执行结果是如下：

x点 x点处的微分值

0.000000 0.000000

0.100000 0.000500

0.200000 0.003000

0.300000 0.009500

0.400000 0.022000

0.500000 0.042504

0.600000 0.073023

0.700000 0.115607

0.800000 0.172408

0.900000 0.245829

1.000000 0.338842

（2）. 只要改动dyfun.m 文件就可以了

function z = fdyfun(x,y)

z = -y + x + 1;

然后取N=20，在MATLAB 命令窗口调用euler.m 文件；

>>[x, y]=euler(@dyfun,[1,1.6],2,20);

执行结果如下：

x 点 x 点处的微分值

1.000000 2.000000

1.030000 2.000450

1.060000 2.001773

1.090000 2.003944

1.120000 2.006937

1.150000 2.010728

1.180000 2.015293

1.210000 2.020610

1.240000 2.026657

1.270000 2.033411

1.300000 2.040852

1.330000 2.048960

1.360000 2.057715

1.390000 2.067097

1.420000 2.077089

1.450000 2.087672

1.480000 2.098829

1.510000 2.110543

1.540000 2.122797

1.570000 2.135575

1.600000 2.148862

2．分别用二分法、不动点迭代法、牛顿法求下列方程的根；

(1)

(2)x 3-x -1=0, x ∈[1,2],ε=10-10（用二分法） x 5-4x -2=0的根，x 0=1.5 （用不动点迭代）

（3）. 用Newton 法求x 3-x -1=0, x ∈[1,2],ε=10-10

解：（1）建立一个m 文件erfen.m, （见附件文件夹）

%二分法

function x=erfen(fun,a,b,ep)

%用途：用二分法求非线性方程f(x)=0有根区间[a,b]中的一个根;

%===============输入=============================%

%格式：x=erfen(fun,a,b,ep)

% fun为函数表达式，

% a, b为区间左右端点，

% ep为精度，

%================================================%

%===============输出=============================%

% x返回近似根

%================================================%

if nargin

error('输入至少有被根函数和根所在区间的上界和下界');

elseif nargin == 3

ep = 1e-15;

end

x=(a+b)/2.0;

i=1;

k=0;

m=abs(subs(fun,x))>ep;

n = (b-a>ep);

while m|n

f1 = subs(fun,a);

f2 = subs(fun,x);

if f1*f2

b=x;

else

a=x;

end

x=(a+b)/2.0;

k=k+1;

X(i)=x;

M(i) = k;

i=i+1;

m=abs(subs(fun,x))>ep;

n = (b-a>ep);

end

disp('迭代次数 对应次数得到的 X 的值');

for i=1:length(X)

fprintf(' %4d %-12f\n',M(i),X(i));

end

disp(['一共迭代的次数是：',num2str(k)]);

disp('最终的近似解是：');

然后建立一个m 文件fun.m （见附件文件夹）

function y = fun(x)

syms x

y = x^3-x-1;

在MATLAB 命令窗口调用erfen.m 文件

>> x=erfen(fun,1,2,1e-10)

执行结果如下：

迭代次数 对应次数得到的 X 的值

1 1.250000

2 1.375000

3 1.312500

4 1.343750

5 1.328125

6 1.320313

7 1.324219

8 1.326172

9 1.325195

10 1.324707

11 1.324951

12 1.324829

13 1.324768

14 1.324738

15 1.324722

16 1.324715

17 1.324718

18 1.324717

19 1.324718

20 1.324718

21 1.324718

22 1.324718

23 1.324718

24 1.324718

25 1.324718

26 1.324718

27 1.324718

28 1.324718

29 1.324718

30 1.324718

31 1.324718

32 1.324718

33 1.324718

34 1.324718

一共迭代的次数是：34

最终的近似解是：

x =

1.3247

（2）把x 5-2或ϕ(x ) =来进行迭代。 x -4x -2=0变为ϕ(x ) =

45

建立m 文件maiter.m(见附件文件夹)

%不动点迭代

function x=maiter(phi,x0,ep,N)

% 用途：用简单迭代法求非线性方程f(x)=0有根区间[a,b]中的一个根 % 格式：x= maiter(phi,x0,ep,N)

% fun为phi(x)的函数句柄,

% x0为初值,

% ep为精度(默认1e-4),

% N为最大迭代次数(默认500),

% x返回近似根

if nargin

if nargin

k=0;

while k

x=feval(phi,x0);

if abs(x-x0)

break;

end

x0=x;k=k+1;

end

if k==N, warning('已达迭代次数上限'); end

disp(['k=',num2str(k)])

再建立一个m 文件fun1.m(见附件文件夹)

function y = fun1(x)

y = (x^5-2)/4;

%y = (4*x+2)^(1/5);

在MATLAB 命令窗口调用maiter.m 文件,N 取10 x 5-2这里用ϕ(x ) =迭代 4

>> x=maiter(@fun1,1.5,10)

k=0

x =

1.3984

这里用ϕ(x ) =

>> x=maiter(@fun1,1.5,10)

k=0

x =

1.5157

（3）建立一个m 文件Newton.m(见附件文件夹)

%Newton迭代法

function x=newton(fun,dfun,x0,ep,N)

%用途：用牛顿法求解非线性方程f(x)=0

%格式：x=newton(fun,dfun,x0,ep,N) fun和dfun 分别为表示f(x)的 %函数句柄, x0为迭代初值, ep为精度(默认1e-4), N为最大迭代 %次数(默认为500), x返回近似根

if nargin

if nargin

k=0;

while k

x=x0-subs(fun,x0)/subs(dfun,x0);

if abs(x-x0)

break;

end

x0=x; k=k+1;

end

if k==N, warning('已达迭代次数上限'); end

disp(['k=',num2str(k)])

再建立一个m 文件dfun.m

function y1 = dfun(x)

syms x

y1 = diff(fun,x);

在MATLAB 命令窗口调用得

>> x=newton(fun,dfun,1,1e-10,20) k=5

x =

1.3247


相关文章

  • 方程求根与解常微分方程

    • 第6章 方程求根与解常微分方程 6.1实验目的 了解微分方程的通解.特解和近似解的概念.熟悉方程求根和常微分方程解的概念,熟悉Mathematica 软件的方程求根和求常微分方程解的命令,掌握用数学软件处理方程求根和常微分方程解的有关问题． ...查看


  • 数值计算基础习题集

    • <数值计算基础>习题集 第1章 引论 1.已知,求近似值的有效数字位数.绝对误差限和相对误差限. 2.下列各数均为四舍五入得到,指出它们各具有几位有效数字及绝对误差限和相对误差限: (1) 6000 (2)7000.00 (3) ...查看


  • 数值计算方法期末考试题

    • y = 1. 已知函数 1 1+x 2的一组数据: 线性插值函数,并计算 求分段 f (1.5) 的近似值. 计算题1. 答案 1 1 ⎰01+x 4. 写出梯形公式和辛卜生公式,并用来分别计算积分. 计算题 4. 答案 四.证明题(本题1 ...查看


  • 算法大全第15章常微分方程的解法

    • 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并加以检验.如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得 ...查看


  • 多种微分方程数值计算方法分析

    • 2010年8月 第4期 城 市 勘 测 Urban Geotechnical Investigation& Surveying Aug. 2010 No. 4 文章编号:1672-8262(2010)04-117-03中图分类号:O ...查看


  • Matlab中常微分方程数值解法

    • Matlab中常微分方程数值解法讲解 作者:dynamic 时间:2008.12.10 版权:All Rights Reserved By www.matlabsky.cn ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ ...查看


  • 实验四 常微分方程初值问题数值解法

    • 实 验 报 告 课程名称 数值分析 实验项目 常微分方程问题初值问题数值解法 一. 实验目的 1.理解如何在计算机上实现用Euler 法.改进Euler 法.Runge -Kutta 算法求一阶常微分方程初值问题 ⎧y '(x ) =f ( ...查看


  • 常微分方程初值问题答案

    • ⎧dy =-y ⎪ 1.(10分)对常微分方程初值问题⎨dx ⎪y (0)=1⎩ (0≤x ≤1) 取步长h =0.1, 分别用改进的Euler 法和标准的四阶Runge-Kutta 法作数值计算, 写出公式和简要推导过程,并把结果填入表内 ...查看


  • 常微分方程常用数值解法

    • 第一章 绪论 1.1 引言 常微分方程是现代数学的一个重要分支,是人们解决各种实际问题的有效工具.微分方程的理论和方法从17世纪末开始发展起来,很快成了研究自然现象的强有力工具,在17到18世纪,在力学.天文.科学技术.物理中,就已借助微分 ...查看


热门内容