欧拉法的matlab程序

欧拉法的matlab程序

f unction oulafa(f,a,b,x0,y0,n,g)

f=inline(f);

g=inline(g);

h=(b-a)/n;

y(1)=y0;

x=a:h:b;

disp('k x(k) y(k) y(x) e ') for i=1:n

y(i+1)=y(i)+h*f(x(i),y(i));

gg(i+1)=g(x(i+1));

e(i+1)=abs(g(i+1)-y(i+1));

fprintf('%d %f %f %f %f\n',i,x(i+1),y(i+1),gg(i+1),e(i+1))

end

matlab实验十七__牛顿迭代法(可打印修改)

实验十七牛顿迭代法 【实验目的】 1.了解牛顿迭代法的基本概念。 2.了解牛顿迭代法的收敛性和收敛速度。 3.学习、掌握MATLAB软件的有关命令。 【实验内容】 用牛顿迭代法求方程的近似根，误差不超过。 3210 ++-=3 10- x x x 【实验准备】 1.牛顿迭代法原理 2.牛顿迭代法的几何解析 3.牛顿迭代法的收敛性 4.牛顿迭代法的收敛速度 5.迭代过程的加速 6.迭代的MATLAB命令 MATLAB中主要用for,while等控制流命令实现迭代。 【实验重点】 1.牛顿迭代法的算法实现 2.牛顿迭代法收敛性和收敛速度 【实验难点】 1.牛顿迭代法收敛性和收敛速度 【实验方法与步骤】 练习1用牛顿迭代法求方程在x=0.5附近的近似 3210 ++-= x x x

根，误差不超过。 310-牛顿迭代法的迭代函数为 322()1()()321 f x x x x g x x x f x x x ++-=-=-'++相应的MATLAB 代码为 >>clear; >>x=0.5; >>for i=1:3 >>x=x-(x^3+x^2+x-1)/(3*x^2+2*x+1) >>end 可算的迭代数列的前3项0.5455，0.5437，0.5437。经三次迭代，就大大超过了精度要求。 练习2 用牛顿迭代法求方程的近似正实根，由此建2(0)x a a =>立一种求平方根的计算方法。 由计算可知，迭代格式为，在实验12的练习4中1()()2a g x x x =+已经进行了讨论。 【练习与思考】 1.用牛顿迭代法求方程的近似根。 ln 1x x =2.为求出方程的根，在区间[1,2]内使用迭代函数进行310x x --=迭代，纪录迭代数据，问迭代是否收敛？对迭代进行加速，对比加速前的数据，比较加速效果。 3.使用在不动点的泰勒公式，证明牛顿迭代法收敛原理。*x

matlab实现牛顿迭代法求解非线性方程组教学文稿

matlab实现牛顿迭代法求解非线性方程组 已知非线性方程组如下 3*x1-cos(x2*x3)-1/2=0 x1^2-81*(x2+0.1)^2+sin(x3)+1.06=0 exp(-x1*x2)+20*x3+(10*pi-3)/3=0 求解要求精度达到0.00001 ———————————————————————————————— 首先建立函数fun 储存方程组编程如下将fun.m保存到工作路径中: function f=fun(x); %定义非线性方程组如下 %变量x1 x2 x3 %函数f1 f2 f3 syms x1 x2 x3 f1=3*x1-cos(x2*x3)-1/2; f2=x1^2-81*(x2+0.1)^2+sin(x3)+1.06; f3=exp(-x1*x2)+20*x3+(10*pi-3)/3; f=[f1 f2 f3]; ———————————————————————————————— 建立函数dfun 用来求方程组的雅克比矩阵将dfun.m保存到工作路径中: function df=dfun(x); %用来求解方程组的雅克比矩阵储存在dfun中 f=fun(x); df=[diff(f,'x1');diff(f,'x2');diff(f,'x3')]; df=conj(df'); ———————————————————————————————— 编程牛顿法求解非线性方程组将newton.m保存到工作路径中: function x=newton(x0,eps,N); con=0; %其中x0为迭代初值eps为精度要求N为最大迭代步数con用来记录结果是否收敛for i=1:N; f=subs(fun(x0),{'x1' 'x2' 'x3'},{x0(1) x0(2) x0(3)}); df=subs(dfun(x0),{'x1' 'x2' 'x3'},{x0(1) x0(2) x0(3)}); x=x0-f/df; for j=1: length(x0); il(i,j)=x(j); end if norm(x-x0)

MATLAB代码 解线性方程组的迭代法

解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;

if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵

牛顿迭代法解元方程组以及误差分析matlab实现

.0],;,[0 ),()(),()(),(0),()(),()(),(,.**,0],;,[),()()(),()()(,0),(),(),(])()[(),(),(),(),(),(])()[(),(),(2,),(])()[(21),(])()[(),(),()(2 )(''))((')()(: 1n 1n 110101010100000000000000000000000000200000000000 00 000fg g f y y g f g f g f fg x x g g f f y x g y y y x g x x y x g y x f y y y x f x x y x f y x y x y x g f g f fg g f y y g f g f g f fg x x g f g f fg g f y y g f g f g f fg x x g g f f y x g y x g y y y x g x x y x f y x f y y y x f x x y x g y x f y x g y y y x x x y x g y x g y x f y x g y x f y y y x x x y x f y x f y x y x f y y y x x x y x f y y y x x x y x f y x f x x f x x x f x f x f x x n n x y y x y y y x y x n n y n n n x n n n n n y n n n x n n n n n x y y x x x x y y x y y x y y x x x x y y x y y y x y x y x y x y y x x y y x x y x y y x x ，则其解可记为： 的行列式不为若系数矩阵： 附近的线性化方程组为在一元方程牛顿迭代法，类似 ，的新近似值于是就得到了根，则可得解： 的行列式不为若系数矩阵），（），（ ），（），（ 则两式构成方程组： 令可得： 构成二元方程组，同样与若另有一方程： 阶小项，得到线性方程忽略在方程根附近取值时，当二元函数的展开为： 开类似一元函数的泰勒展?????+-+=-+-+=?????=-+-+=-+-+??? ????-+-+=-+-+=????????-+-=--+-=-?????-=-+--=-+-==??-+??-+=??-+??-+=??-+??-+??-+??-+=-+ -+=++========η ξξ

实验五 欧拉法Matlab实验报告

北京理工大学珠海学院实验报告 ZHUHAI CAMPAUS OF BEIJING INSTITUTE OF TECHNOLOGY 班级2012电气2班学号120109021010姓名陈冲指导教师张凯成绩 实验题目（实验五）欧拉法实验地点及时间JD501 2014/1/2（6-7节） 一、实验目的 1.掌握用程序语言来编辑函数。 2.学会用MATLAB编写Euler.m以及TranEuler.m函数。 二、实验环境 Matlab软件 三、实验内容 1、以书中第124页题目11为例编辑程序来实现计算结果。 2、使用MATLAB进行编写： 第一步：编写Euler.m函数，代码如下 编写TranEuler.m函数，代码如下 第二步：利用上述函数编辑命令：（可见实验结果中的截图）

在此之前先建立一个名为f.m 的M 文件，代码如下 function z=f(x); z=8-3y; 再编辑代码： 得到了欧拉法的结果：y （0.4）=2.47838030901267 编辑另一段命令： 得到改进欧拉法的结果：y （0.4）=2.46543714659780 在此基础上，我还编辑龙格库达的命令窗口代码，如下： 四、实验题目 用欧拉法和改进欧拉法求解初值问题'83,(0)2y y y =-=，试取步长0.2h =计算(0.4)y 的近似值。 五、实验结果

六、总结 通过这次实验我掌握了将得到的解进一步精确，而且要学会比较这几种方法的精确性，显然，四阶龙格库达比改进欧拉发精确，改进欧拉发比欧拉法精确。 实验难度不大，要比较n的取值不同，产生的影响不同。

MATLAB程序(牛顿法及线形方程组)

MATLAB 程序 1、图示牛顿迭代法（M 文件）文件名：newt_g function x = new_g(f_name,x0,xmin,xmax,n_points) clf,hold off % newton_method with graphic illustration del_x = 0.001; wid_x = xmax - xmin; dx = (xmax - xmin)/n_points; xp = xmin:dx:xmax; yp = feval(f_name,xp); plot(xp,yp);xlabel('x');ylabel('f(x)'); title('newton iteration'),hold on ymin = min(yp); ymax = max(yp); wid_y = ymax-ymin; yp = 0. * xp; plot(xp,yp) x = x0; xb = x+999; n=0; while abs(x-xb) > 0.000001 if n > 300 break; end y=feval(f_name,x); plot([x,x],[y,0]);plot(x,0,'o') fprintf(' n = % 3.0f, x = % 12.5e, y = % 12.5e \ n', n, x, y); xsc = (x-xmin)/wid_x; if n < 4, text(x,wid_y/20,[num2str(n)]), end y_driv = (feval(f_name,x + del_x) - y)/del_x; xb = x; x = xb - y/y_driv; n = n+1; plot([xb,x],[y,0]) end plot([x x],[0.05 * wid_y 0.2 * wid_y]) text( x, 0.2 * wid_y, 'final solution') plot([ x ( x - wid_x * 0.004)], [0.01 * wid_y 0.09 * wid_y]) plot([ x ( x + wid_x * 0.004)], [0.01 * wid_y 0.09 * wid_y]) 传热问题 假设一个火炉是用厚度为0.05m 的砖单层砌成的。炉内壁温度为T 0=625K, 外壁温度为T 1（未知）。由于对流和辐射造成了外壁的热量损失，温度T 1由下式决定： 44111()()()()0f k f T T T T T h T T x εσ∞=-+-+-=? 其中： k ：炉壁的热传导系数，1.2W/mK ε: 发射率，0.8 T 0：内壁温度，625K T 1：外壁温度（未知），K T ∞：环境温度，298K T f ：空气温度，298K H ：热交换系数，20W/m 2K

欧拉法matlab程序

法 function [x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler(dyfun,[0,1],1,;[x,y] ans = 2.隐式Euler法 function [x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if k>K error('迭代发散'); end end >> dyfun=inline('y-2*x/y');

[x,y]=naeulerb(dyfun,[0,1],1,;[x,y] ans = 3.改进Euler法 function [x,y]=naeuler2(dyfun,xspan,y0,h) 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'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler2(dyfun,[0,1],1,;[x,y] ans =

lu分解法、列主元高斯法、jacobi迭代法、gaussseidel法的原理及matlab程序

一、实验目的及题目 1.1 实验目的： （1）学会用高斯列主元消去法，LU 分解法，Jacobi 迭代法和Gauss-Seidel 迭代法解线性方程组。 （2）学会用Matlab 编写各种方法求解线性方程组的程序。 1.2 实验题目： 1. 用列主元消去法解方程组： 1241234 123412343421233234x x x x x x x x x x x x x x x ++=??+-+=??--+=-??-++-=? 2. 用LU 分解法解方程组,Ax b =其中 4824012242412120620266216A --?? ?- ?= ? ?-??，4422b ?? ? ?= ?- ?-?? 3. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组： 123234 1231234102118311210631125x x x x x x x x x x x x x -+=-??-+=-??-+=??-+-+ =? 二、实验原理、程序框图、程序代码等 2.1实验原理 2.1.1高斯列主元消去法的原理 Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式： 1111221122222n n n n nn n n b x b x b x g b x b x g b x g +++=??++=????= ? 这个过程就是消元，然后再回代就好了。具体过程如下： 对于1,2, ,1k n =-，若() 0,k kk a ≠依次计算

()() (1)()()(1)()()/,,1, ,k k ik ik kk k k k ij ij ik kj k k k i i ik k m a a a a m a b b m b i j k n ++==-=-=+ 然后将其回代得到： ()() ()()()1/()/,1,2,,1 n n n n nn n k k k k k kj j kk j k x b a x b a x a k n n =+?=??=-=--? ? ∑ 以上是高斯消去。 但是高斯消去法在消元的过程中有可能会出现() 0k kk a =的情况，这时消元就无法进行了，即使主元数() 0,k kk a ≠但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。 2.1.2直接三角分解法（LU 分解）的原理 先将矩阵A 直接分解为A LU =则求解方程组的问题就等价于求解两个三角形方程组。 直接利用矩阵乘法，得到矩阵的三角分解计算公式为： 1111111 11 1,1,2,,/,2,,,,,1,,,2,3, ()/,1,2, ,i i i i k kj kj km mj m k ik ik im mk kk m u a i n l a u i n u a l u j k k n k n l a l u u i k k n k n -=-===?? ==?? =-=+??=??=-=++≠?? ∑∑且 由上面的式子得到矩阵A 的LU 分解后，求解Ux=y 的计算公式为 11 111,2,3,/()/,1,2, ,1 i i i ij j j n n nn n i i ij j ii j i y b y b l y i n x y u x y u x u i n n -==+=??? =-=?? =??? =-=--?? ∑∑ 以上为LU 分解法。

matlab 欧拉算法 附截图

设系统方程为：y t y y /2)1(-=,1)0(=y ,用改进欧拉法求解各离散点y 的数值解，步长 10,1.0≤≤=t h ，解析解为t y 21+= 。 解：改进欧拉法 ),(1n n n p n y t hf y y +=+ )],(),([5.0111p n n n n n c n y t f y t f h y y +++++= 已知 n n n n n y t y y t f /2),(-= n n n n n n n p n y ht y h y t y h y y /2)1()/2(1-+=-+=+ 1 111111/5.0/)5.01()]/2()/2[(5.0+++++++-+-+=-+-+=n n n n n n n n n n n n n c n y ht hy y ht y h y t y y t y h y y 程序： h=0.1; t=0:h:1; N=length(t); y=ones(1,N); ey=ones(1,N); zy=ones(1,N); for k=1:N-1 y(1,k+1)=(1+h)*y(1,k)-(2*h*(k-1)/(N-1))./y(1,k);%预估公式 ey(1,k+1)=(1+h)*ey(1,k)-(2*h*(k-1)/(N-1))./ey(1,k);%欧拉公式 y(1,k+1)=(1+0.5*h)*y(1,k)-(h*(k-1)/(N-1))./y(1,k)+0.5*h*y(1,k+1)-(h*k/(N-1))./y(1,k+1);%改进欧拉 zy(1,k+1)=(1+2*k/(N-1)).^0.5;%解析解 end plot(t,zy,'-xk',t,y,':ob',t,ey,'-.*r','linewidth',1.0); xlabel('t'); ylabel('y'); 截图：

数值分析的MATLAB程序

列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);

end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')

MATLAB样例之雅克比迭代法

要求： 下面分别使用雅克比迭代法和高斯-赛德尔迭代法求一个方程组的近似解用的线性方程组是按实验要求给的： 7*x1+x2+2*x3=10 x1+8*x2+2*x3=8 2*x1+2*x2+9*x3=6 雅克比迭代法的matlab代码：（老师写的） A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(any(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); while 1 x1=B*x0+f K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf)N fprintf('迭代超限') end x0=x1; end 高斯-赛德尔迭代法matlab代码:（自己改的）

A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(all(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); x00=x0; while 1 x11=B*x0+f; x00(1,1)=x11(1,1); x12=B*x00+f; x00(2,1)=x12(2,1); x13=B*x00+f; x00(3,1)=x13(3,1); x1=x00 K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf)N fprintf('迭代超限') end x0=x1; end

Euler法解微分方程-Matlab程序

%主程序main.m-----OK! clear; T=0.0; Y=zeros(3,1); Y(1)=1.0;Y(2)=1.0;Y(3)=1.0; H1=0.05;M=3;EPS=1.0e-05;%EPS精度要求M方程个数H1拟定的输出步长 for i=1:10 [X,Y]=euler1(T,H1,Y,M,EPS) T=T+H1; End %变步长euler方法 function [X,Y1]=euler1(T,H1,Y,M,EPS) %M-方程个数,EPS-精度,Y0-右端初值，T-自变量前一点值，H-步长 N=1;P=1+EPS;X=T;G=zeros(M,1); H=H1;%H-在程序中要改变的步长H1－主程序中确定的输出步长 for i=1:M C(i)=Y(i); end K1=zeros(M,1);K2=zeros(M,1);K3=zeros(M,1);K4=zeros(M,1); while P>=EPS %变步长积分一步（H1） for i=1:M G(i)=Y(i); Y(i)=C(i); end DT=H/N; T=X; %－－变步长积分过程 for j=1:N K1=F(Y); K2=F(Y+H/2*K1'); K3=F(Y+H/2*K2'); K4=F(Y+H*K3'); for i=1:M Y(i)=Y(i)+H/6*(K1(i)+2*K2(i)+2*K3(i)+K4(i)); T=T+DT; end end %－－－－－－－－－－－－－－－－－－－－－ P=0.0; for i=1:M Q=abs(Y(i)-G(i)); if Q>P

P=Q; end end H=H/2.0; N=N+N; end T=X; X=T+H1; Y1=Y; %右端函数值function D=F(y) D(1)=y(2); D(2)=-1*y(1); D(3)=y(3);

非线性方程组求解的牛顿迭代法用MATLAB实现

1. 二元函数的newton 迭代法理论分析 设),(y x f z =在点),(00y x 的某一邻域内连续且有直到2阶的连续偏导数，),(00h y h x ++为该邻域内任意一点，则有 ?? ? ????? +??+≈++==00) ,(),(),(),(0000y y x x y x f y k y x f x h y x f k y h x f 其中 0x x h -=，0y -=y k 于是方程0),(=y x f 可近似表示为 0) ,(),(),(k =?? ? ????? +??+==k k y y x x k y x f y k y x f x h y x f 即 0),()(),()(),(y k =-+-+k k k k k x k k y x f y y y x f x x y x f 同理，设y)g(x,z =在点),(00y x 的某一邻域内连续且有直到2阶的连续偏导数，),(00h y h x ++为该邻域内任意一点，亦有 ?? ?????? +??+≈++==00),(),(),(),(0000y y x x y x g y k y x g x h y x g k y h x g 其中0x x h -=，0y -=y k 于是方程0),(g =y x 可近似表示为 0) ,(),(),(k =?? ? ????? +??+==k k y y x x k y x g y k y x g x h y x g 即 0),(g )(),()(),(y k =-+-+k k k k k x k k y x y y y x g x x y x g 于是得到方程组 ? ??=-+-+=-+-+0),(g )(),()(),(0),()(),()(),(y k y k k k k k k x k k k k k k k x k k y x y y y x g x x y x g y x f y y y x f x x y x f

MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程

姓名：樊元君学号：2012200902 日期：2012.11.06 一、实验目的 掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。掌握使用MATLAB程序求解常微分方程问题的方法。 二、实验内容 1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。 实验中以下列数据验证程序的正确性。 求，步长h=0.25。 2、实验注意事项 的精确解为，通过调整步长，观察结果的精度的变化

三、程序流程图： ●改进欧拉格式流程图：

●四阶龙格库塔流程图： 四、源程序： ●改进后欧拉格式程序源代码：

format long h=input('h='); x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); X=input('X=');Y=input('Y='); n=round((Y-X)/h); i=1;x1=0;yp=0;yc=0; for i=1:1:n x1=x0+h; yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);% yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);% y1=(yp+yc)/2; x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y); end end ●四阶龙格库塔程序源代码：

2-8牛顿迭代法matlab

实验七 牛顿迭代法 【实验目的】 1．了解牛顿迭代法的基本概念。 2．了解牛顿迭代法的收敛性和收敛速度。 3．学习掌握MATLAB 软件有关的命令。 【实验内容】 用牛顿迭代法求方程0123=-++x x x 的近似根，误差不超过310-。 【实验准备】 1．牛顿迭代法原理 设已知方程0)(=x f 的近似根0x ,则在0x 附近)(x f 可用一阶泰勒多项式))((')()(000x x x f x f x p -+=近似代替.因此, 方程0)(=x f 可近似地表示为0)(=x p .用1x 表示0)(=x p 的根,它与0)(=x f 的根差异不大. 设0)('0≠x f ,由于1x 满足,0))((')(0100=-+x x x f x f 解得 ) (')(0001x f x f x x -= 重复这一过程,得到迭代格式 ) (')(1n n n n x f x f x x -=+ 这就是著名的牛顿迭代公式,它相应的不动点方程为 ) (')()(x f x f x x g -=. 2. 牛顿迭代法的几何解析 在0x 处作曲线的切线，切线方程为))((')(000x x x f x f y -+=。令 0=y ，可得切线与x 轴的交点坐标) (')(0001x f x f x x -=，这就是牛顿法的迭代公式。因此，牛顿法又称“切线法”。

3．牛顿迭代法的收敛性 计算可得2)] ('[)(")()('x f x f x f x g -=,设*x 是0)(=x f 的单根,有0)(',0)(**≠=x f x f ,则 0)]('[)(")()('2**** =-=x f x f x f x g , 故在*x 附近,有1)('>clear; >>x=0.5; >>for i=1:3 >>x=x-(x^3+x^2+x-1)/(3*x^2+2*x+1) >>end 可算得迭代数列的前3项0.5455, 0.5437, 0.5437.近三次迭代,就大大超过了精度要求. 练习２用牛顿迭代法求方程)0(2>=a a x .的近似正实根,由此建立一种求平方根的计算方法. 由计算可知,迭代格式为)(21)(x a x x g += .,在实验12的练习4种已经进行了讨论. 练习3用牛顿迭代法求方程1=x xe 的正根. 牛顿迭代法的迭代函数为

二分法、简单迭代法的matlab代码实现

实验一非线性方程的数值解法(一) 信息与计算科学金融崔振威201002034031一、实验目的： 熟悉二分法和简单迭代法的算法实现。 二、实验内容： 教材P40 2.1.5 三、实验要求 1根据实验内容编写二分法和简单迭代法的算法实现 2简单比较分析两种算法的误差 3试构造不同的迭代格式，分析比较其收敛性 (一)、二分法程序： function ef=bisect(fx,xa,xb ,n, delta) % fx是由方程转化的关于x的函数，有fx=0。 % xa解区间上限 % xb解区间下限 % n最多循环步数，防止死循环。 %delta为允许误差 x=xa;fa=eval(fx); x=xb;fb=eval(fx); disp(' [ n xa xb xc fc ]'); for i=1: n xc=(xa+xb)/2;x=xc;fc=eval(fx); X=[i,xa,xb,xc,fc]; disp(X), if fc*fa<0 xb=xc; else xa=xc; end if (xb-xa)

k=0; while abs(x-xO)>eps & k> fplot('[x A5-3*x A3-2*x A2+2]',[-3,3]);grid 得下图： 由上图可得知：方程在［-3,3］区间有根。 （2 ）、二分法输出结果 >> f='xA5-3*xA3-2*xA2+2' f = X A5-3*X A3-2*X A2+2 >> bisect(f,-3,3,20,10A(-12)) 2.0000 - 3.0000 0 -1.5000 0.0313

MATLAB Euler法解常微分方程

Euler 法解常微分方程 Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算h n n +=判断b n ≤是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算),(1n n n n y x hf y y +=+ Step 4 ),(1n n n n y x hf y y +=+ Euler 法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

0.4613 指导教师： 年 月 日 改进Euelr 法解常微分方程 改进Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2 取一个以h 为步长，a ，b 分别为左右端点的矩阵 Step 3 （1）做显性Euler 预测),( 1n n i i y x hf y y +=+ （2）将1+i y 带入,(),([2h 11++++=i i i i i x f y x f y y Step 4计算h n n +=判断b n ≤是否成立，成立返回Step 5 )],(),([2h 111+++++=i i i i i i y x f y x f y y 改进Euler 法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x'

matlab 迭代法[精品]

matlab 迭代法[精品] 1. 矩阵 122,211,，，，，,,,,A,111A,222, 11,,,,,,,,221,,112，，，， 证明:求解以为系数矩阵线性方程组的Jacobi迭代式收敛的，而A1 Gauss-Seidel方法是发散的;求解以为系数矩阵线性方程组的A2实验名称Gauss-Seidel是收敛的，而Jacobi方法是发散的. 2. 矩阵 1aa，，,,Aaa,1 ,,,,aa1，， (a) 参数取什么值时，矩阵是正定的. a (b) 取什么值时，求以为系数矩阵线性方程组的Jacobi迭代式收aa 敛的. 1、根据迭代收敛性的充分必要条件来判断Jacobi迭代式与Gauss-Seide 迭代式的收敛性，迭代收敛性仅与方程组系数矩阵有关，与右端无关;而且不依赖于初值的选取。实验目的 2、根据矩阵的判断定理求得矩阵元素a的取值，同时根据矩阵线性方程组的Jacobi迭代式收敛的充分条件(严格对角占优)来求a得取值。 1、(1)检验线性方程组的Jacobi迭代式的收敛性: function jacobi(A) D=zeros(3); for i=1:3 D(i,i)=A(i,i); 实验内容end (算法、程B=D^(-1)*(D-A); 序、步骤和k=max(abs(eig(B))) 方法) if k<1

'该线性方程组的Jacobi迭代式是收敛的' else k>=1 '该线性方程组的Jacobi迭代式是发散的' end (2)检验线性方程组的Gauss-Seide迭代式的收敛性: function Gauss(A) D=zeros(3); L=zeros(3); U=zeros(3); for i=1:3 D(i,i)=A(i,i); end L(2:3,1)=A(2:3,1); L(3,2)=A(3,2); U(1,2:3)=A(1,2:3); U(2,3)=A(2,3); B=-(D+L)^(-1)*U; k=max(abs(eig(B))) if k<1 '该线性方程组的Gauss-Seidel迭代式是收敛的' else k>=1 '该线性方程组的Gauss-Seidel迭代式是发散的' end 2、(1)参数取什么值时，矩阵是正定的.(矩阵的特征值全为正) a >> syms a >> A=[1 a a;a 1 a;a a 1]; >> eig(A) ans = 2*a+1 1-a