matlab实验 非线性方程(组)求解
- 格式:docx
- 大小:124.84 KB
- 文档页数:11
非线性方程组求解-Matlab-fsolve实例一:①建立文件fun.m:function y=fun(x)y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ...x(2) - 0.5*cos(x(1))+0.3*sin(x(2))];②>>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve'))注:...为续行符m文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。
实例二:①建立文件fun.mfunction F=myfun(x)F=[x(1)-3*x(2)-sin(x(1));2*x(1)+x(2)-cos(x(2))];②然后在命令窗口求解:>> x0=[0;0]; %设定求解初值>> options=optimset('Display','iter'); %设定优化条件>> [x,fv]=fsolve(@myfun,x0,options) %优化求解%MATLAB显示的优化过程Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius0 3 1 2 11 6 0.000423308 0.5 0.0617 12 9 5.17424e-010 0.00751433 4.55e-005 1.253 12 9.99174e-022 1.15212e-005 9.46e-011 1.25 Optimization terminated: first-order optimality is less than options.TolFun.x =0.49660.0067fv =1.0e-010 *0.31610.0018实例三:求下列非线性方程组在(0.5,0.5) 附近的数值解。
1大学数学实验 实验报告 | 2014/4/5一、 实验目的1、学习用Matlab 软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2、通过实例学习用线性代数方程组解决简化问题。
二、 实验内容项目一:种群的繁殖与稳定收获:种群的数量因繁殖而增加,因自然死亡而减少,对于人工饲养的种群(比如家畜)而言,为了保证稳定的收获,各个年龄的种群数量应维持不变。
种群因雌性个体的繁殖而改变,为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n ,当年年龄k 的种群数量记作x k ,繁殖率记作b k (每个雌性个体1年的繁殖的数量),自然存活率记作s k (s k =1−d k ,d k 为1年的死亡率),收获量记作ℎk ,则来年年龄k 的种群数量x ̌k 应该为x ̌k =∑b k n k=1x k , x ̌k+1=s k x k −ℎk , (k=1,2,…,n -1)。
要求各个年龄的种群数量每年维持不变就是要求使得x ̌k =x k , (k=1,2,…,n -1).(1) 如果b k , s k 已知,给定收获量ℎk ,建立求各个年龄的稳定种群数量x k 的模型(用矩阵、向量表示).(2) 设n =5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如要求ℎ1~ℎ5为500,400,200,100,100,求x 1~x 5.(3) 要使ℎ1~ℎ5均为500,如何达到?问题分析:该问题属于简单的种群数量增长模型,在一定的条件(存活率,繁殖率等)下为使各年龄阶段的种群数量保持不变,各个年龄段的种群数量将会满足一定的要求,只要找到种群数量与各个参量之间的关系,建立起种群数量恒定的方程就可以求解出各年龄阶段的种群数量。
模型建立:根据题目中的信息,令x ̌k =x k ,得到方程组如下:{x ̌1=∑b k nk=1x k =x 1x ̌k+1=s k x k −ℎk =x k+1整理得到:{−x 1∑b k nk=1x k =0−x k+1+s k x k =ℎk2 大学数学实验 实验报告 | 2014/4/52写成系数矩阵的形式如下:A =[b 1−1b 2b 3s 1−100s 2−1…b n−1b n0000⋮⋱⋮000000000⋯00−10s n−1−1]令h =[0, ℎ1,ℎ2,ℎ3,…,ℎn−2,ℎn−1]Tx =[x n , x n−1,…,x 1]T则方程组化为矩阵形式:Ax =h ,即为所求模型。
实验六非线性方程求解姓名:李严凯学号:2009011887 班级:化93实验目的1掌握用MATLAB 软件求解非线性方程和方程组的基本用法,并对结果作出初步的分析2 练习用非线性方程组建立实际问题的模型并进行求解实验内容3.问题复述:(1)小张夫妇以按揭方式贷款买了一套价值20万的房子,首付了5万元,每月还款1000元,15年还清。
问贷款利率是多少?(2)某人欲贷款50万元购房,他咨询了两家银行,第一家银行开出的条件是每月还4500元,15年还清;第二家银行开出的条件是每年还45000元,20年还清。
从利率方面看,那家银行较优惠(简单假设年利率=月利率*12)?解答:模型:在实际问题中,购房人在支付首付后,向银行按揭申请贷款,并由银行支付购房人不足以支付的部分房价。
而后购房人将在一定时期内定期定额缴纳月供,直至还清银行的贷款本金及利息。
而银行则每月设定固定的月利率,按照复利的方式计算购房人应需缴纳的房款。
不妨设购房之初,购房人向银行申请了a0的贷款,银行的月利率为x,每月购房人向银行支付的月供为b,则在以后各月(第1、2、3…n个月),购房人所欠银行的贷款为:a1=a0*(1+x)-b;a2=a1*(1+x)-b;a3=a2*(1+x)-b;……an=a(n-1)*(1+x)-b当购房人所欠银行的贷款数an=0时,其贷款还清,还贷过程结束。
到此模型分析完毕,可以进行MATLAB求解。
求解:首先写出函数,表示出第n个月购房人所欠银行贷款:function y = interest( x,n ) %x为月利率n=180; %n为还清贷款的月数a(1)=15*(1+x)-0.1; %第1个月欠款for i=2:na(i)=a(i-1)*(1+x)-0.1; %第n个月欠款endy=a(n);end得到以上表达式后,令其中y=0即可解得月利率x,编程如下:[x,fv,ef,out]=fzero(@interest,0.05)考虑到银行贷款利率的实际情况,初值设为了0.05程序输出结果如下:x = 0.002081163889460fv = -5.287437154777308e-014ef =1out = intervaliterations: 12iterations: 14funcCount: 38algorithm: 'bisection, interpolation'message: 'Zero found in the interval [-0.014, 0.0952548]' 当然,也可以采用fsolve求解,列出:[x,fv,ef,out]=fsolve(@interest,0.05)得到x =0.002081163889460fv =-5.287437154777308e-014ef =1out = iterations: 14funcCount: 30algorithm: 'trust-region dogleg'firstorderopt: 9.769300386433235e-011message: [1x695 char]可以看出,两种解法的结果基本相同,不过fzero的解答时间更快,效果更优。
matlab求解非线性方程组及极值默认分类2010-05-18 15:46:13 阅读1012 评论2 字号:大中小订阅一、概述:求函数零点和极值点:Matlab中三种表示函数的方法: 1. 定义一个m函数文件, 2.使用函数句柄; 3.定义inline函数, 其中第一个要掌握简单函数编写, 二, 三中掌握一个。
函数的'常规'使用有了函数了, 我们怎么用呢, 一种是直接利用函数来计算, 例如: sin(pi), 还有我们提到的mysqr(3)...另一种是函数画图, 例如Plottools中提到的ezplot, ezsurf... 但是这也太小儿科了, 有没有想过定义函数后, 利用它来: 求解零点(即解f(x)=0方程), 最优化(求最值/极值点), 求定积分, 常微分方程求解等. 当然这里由于篇幅有限(空间快满了)以及这个只是'基础教程'的缘故, 只提及一些皮毛知识, 掌握这些后, 如果需要你可以进一步学习.解f(x)=0已知函数求解函数值=0所表示的方程, Matlab中有两个函数可以做到, fzero和fsolve前者只能解一元方程, 后者可以解多元方程组, 不过基本使用形式上差不多:解=fzero(函数, 初值, options)解=fsolve(函数, 初值, options)关于解: fzero给出的是x单值的解, fsolve给出的是解x可能处于的区间, 当然, 这个区间很窄.关于'函数', 还记得前面提到的三种表示方法吧, 在这里都可以用, 记住就是: 如果直接使用函数名, 要用单引号将它括起来, 而函数句柄, inline函数可以直接使用.关于'初值': 电脑比较笨, 它寻找解的办法是尝试不同地x值, 摸索解在哪里, 所以我们一开始就要给它指明从哪里开始下手, 初值这里, 可以只给它一个值, 让它在这个值附近找解, 也可以给它一个区间(区间用[下限,上限]这种方式表示), 它会在这个区间内找解.fzero的一些局限, 如果你给定的初值是区间, 而恰好函数在区间端点处同号, fzero会出错, 而如果你只给一个初值, fezro又有可能'走错方向', 例如给初值2让它解mysqr这个函数方程就出错了, FT!寻找函数极值/最值Matlab中也有两个函数可以做到, 是: fminbnd: 寻找一元函数极小值; fminsearch: 寻找多元函数极小值(当然一元也行). 别问我怎么没有找极大值的Matlab函数, 你把原函数取负数, 寻找它的极小值不就行了. 相关语法:x=fminbnd(函数, 区间起始值, 区间终止值)x=fminsearch(函数, 自变量初值)相关说明: fminbnd中指定要查找极小值的自变量区间, 好像不指定也行, 不过那样的话, 如果函数有多个极小值就可能比较难以预料结果了.fminsearch中要给定一个初值, 这个初值可以是自变量向量(将自变量依次排在一起组成向量)的初值, 也可以是表示向量初值区间的一个矩阵.函数: 那三种形式都适用, 但是记住, 直接使用函数名称需要加单引号!cite from:/qq529312840/blog/item/3687e4c7e7e2d6d9d0006049.html二、实例+讲解(1)非线性方程数值求解:1 单变量非线性方程求解在MATLAB中提供了一个fzero函数,可以用来求单变量非线性方程的根。
非线性方程组求解1.mulStablePoint用不动点迭代法求非线性方程组的一个根function [r,n]=mulStablePoint(F,x0,eps)%非线性方程组:f%初始解:a%解的精度:eps%求得的一组解:r%迭代步数:nif nargin==2eps=1.0e-6;endx0 = transpose(x0);n=1;tol=1;while tol>epsr= subs(F,findsym(F),x0); %迭代公式tol=norm(r-x0); %注意矩阵的误差求法,norm为矩阵的欧几里德范数n=n+1;x0=r;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endend2.mulNewton用牛顿法法求非线性方程组的一个根function [r,n]=mulNewton(F,x0,eps)if nargin==2eps=1.0e-4;endx0 = transpose(x0);Fx = subs(F,findsym(F),x0);var = findsym(F);dF = Jacobian(F,var);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx;n=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx; %核心迭代公式tol=norm(r-x0);n=n+1;if(n>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endend3.mulDiscNewton用离散牛顿法法求非线性方程组的一个根function [r,m]=mulDiscNewton(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=transpose(x0)-inv(J)*fx;m=1;tol=1;while tol>epsxs=r;fx = subs(F,findsym(F),xs);J = zeros(n,n);for i=1:nx1 = xs;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=xs-inv(J)*fx; %核心迭代公式tol=norm(r-xs);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;4.mulMix用牛顿-雅可比迭代法求非线性方程组的一个根function [r,m]=mulMix(F,x0,h,l,eps)if nargin==4eps=1.0e-4;endn = length(x0);J = zeros(n,n);Fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));C =D - J;inD = inv(D);H = inD*C;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = Hm*inD*Fx;r = transpose(x0)-dr; m=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));C =D - J;inD = inv(D);H = inD*C;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = Hm*inD*Fx;r = x0-dr; %核心迭代公式tol=norm(r-x0);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endend5.mulNewtonSOR用牛顿-SOR迭代法求非线性方程组的一个根function [r,m]=mulNewtonSOR(F,x0,w,h,l,eps)if nargin==5eps=1.0e-4;endn = length(x0);J = zeros(n,n);Fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));L = -tril(J-D);U = -triu(J-D);inD = inv(D-w*L);H = inD*(D - w*D+w*L);;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = w*Hm*inD*Fx;r = transpose(x0)-dr;m=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));L = -tril(J-D);U = -triu(J-D);inD = inv(D-w*L);H = inD*(D - w*D+w*L);;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = w*Hm*inD*Fx;r = x0-dr; %核心迭代公式tol=norm(r-x0);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endend6.mulDNewton用牛顿下山法求非线性方程组的一个根function [r,m]=mulDNewton(F,x0,eps)%非线性方程组:F%初始解:x0%解的精度:eps%求得的一组解:r%迭代步数:nif nargin==2eps=1.0e-4;endx0 = transpose(x0);dF = Jacobian(F);m=1;tol=1;while tol>epsttol=1;w=1;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);F1=norm(Fx);while ttol>=0 %下面的循环是选取下山因子w的过程r=x0-w*inv(dFx)*Fx; %核心的迭代公式Fr = subs(F,findsym(F),r);ttol=norm(Fr)-F1;w=w/2;endtol=norm(r-x0);m=m+1;x0=r;if(m>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endend7.mulGXF1用两点割线法的第一种形式求非线性方程组的一个根function [r,m]=mulGXF1(F,x0,x1,eps)format long;if nargin==3eps=1.0e-4;endx0 = transpose(x0);x1 = transpose(x1);n = length(x0);fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);for i=1:nxt = x1;xt(i) = x0(i);J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);endr=x1-inv(J)*fx1;m=1;tol=1;while tol>epsx0 = x1;x1 = r;fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);for i=1:nxt = x1;xt(i) = x0(i);J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);endr=x1-inv(J)*fx1;tol=norm(r-x1);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;8.mulGXF2用两点割线法的第二种形式求非线性方程组的一个根function [r,m]=mulGXF2(F,x0,x1,eps)format long;if nargin==3eps=1.0e-4;endx0 = transpose(x0);x1 = transpose(x1);n = length(x0);fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);xt = x1;xt(1) = x0(1);J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);for i=2:nxt = x1;xt(1:i) = x0(1:i);xt_m = x1;xt_m(1:i-1) = x0(1:i-1);J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);endr=x1-inv(J)*fx1;m=1;tol=1;while tol>epsx0 = x1;x1 = r;fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);xt = x1;xt(1) = x0(1);J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);for i=2:nxt = x1;xt(1:i) = x0(1:i);xt_m = x1;xt_m(1:i-1) = x0(1:i-1);J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);endr=x1-inv(J)*fx1;tol=norm(r-x1);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;9.mulVNewton用拟牛顿法求非线性方程组的一组解function [r,m]=mulVNewton(F,x0,A,eps)%方程组:F%方程组的初始解:x0% 初始A矩阵:A%解的精度:eps%求得的一组解:r%迭代步数:mif nargin==2A=eye(length(x0)); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendx0 = transpose(x0);Fx = subs(F, findsym(F),x0);r=x0-A\Fx;m=1;tol=1;while tol>epsx0=r;Fx = subs(F, findsym(F),x0);r=x0-A\Fx;y=r-x0;Fr = subs(F, findsym(F),r);z= Fr-Fx;A1=A+(z-A*y)*transpose(y)/norm(y); %调整A A=A1;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end10.mulRank1用对称秩1算法求非线性方程组的一个根function [r,n]=mulRank1(F,x0,A,eps)if nargin==2l = length(x0);A=eye(l); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-inv(A)*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-inv(A)*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;A1=A+ fr *transpose(fr)/(transpose(fr)*y); %调整A A=A1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end11.mulDFP用D-F-P算法求非线性方程组的一组解function [r,n]=mulDFP(F,x0,A,eps)if nargin==2l = length(x0);B=eye(l); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-B*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-B*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;B1=B+ y*y'/(y'*z)-B*z*z'*B/(z'*B*z); %调整AB=B1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end12.mulBFS用B-F-S算法求非线性方程组的一个根function [r,n]=mulBFS(F,x0,B,eps)if nargin==2l = length(x0);B=eye(l); %B取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-B*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-B*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;u = 1 + z'*B*z/(y'*z);B1= B+ (u*y*y'-B*z*y'-y*z'*B)/(y'*z); %调整B B=B1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end13.mulNumYT用数值延拓法求非线性方程组的一组解function [r,m]=mulNumYT(F,x0,h,N,eps)format long;if nargin==4eps=1.0e-8;endn = length(x0);fx0 = subs(F,findsym(F),x0);x0 = transpose(x0);J = zeros(n,n);for k=0:N-1fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endinJ = inv(J);r=x0-inJ*(fx-(1-k/N)*fx0);x0 = r;endm=1;tol=1;while tol>epsxs=r;fx = subs(F,findsym(F),xs);J = zeros(n,n);for i=1:nx1 = xs;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=xs-inv(J)*fx; %核心迭代公式tol=norm(r-xs);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;14.DiffParam1用参数微分法中的欧拉法求非线性方程组的一组解function r=DiffParam1(F,x0,h,N)%非线性方程组:f%初始解:x0%数值微分增量步大小:h%雅可比迭代参量:l%解的精度:eps%求得的一组解:r%迭代步数:nx0 = transpose(x0);n = length(x0);ht = 1/N;Fx0 = subs(F,findsym(F),x0);for k=1:NFx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endinJ = inv(J);r = x0 - ht*inJ*Fx0;x0 = r;end15.DiffParam2用参数微分法中的中点积分法求非线性方程组的一组解function r=DiffParam2(F,x0,h,N)%非线性方程组:f%初始解:x0%数值微分增量步大小:h%雅可比迭代参量:l%解的精度:eps%求得的一组解:r%迭代步数:nx0 = transpose(x0);n = length(x0);ht = 1/N;Fx0 = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nxt = x0;xt(i) = xt(i)+h(i);J(:,i) = (subs(F,findsym(F),xt)-Fx0)/h(i);endinJ = inv(J);x1 = x0 - ht*inJ*Fx0;for k=1:Nx2 = x1 + (x1-x0)/2;Fx2 = subs(F,findsym(F),x2);J = zeros(n,n);for i=1:nxt = x2;xt(i) = xt(i)+h(i);J(:,i) = (subs(F,findsym(F),xt)-Fx2)/h(i);endinJ = inv(J);r = x1 - ht*inJ*Fx0;x0 = x1;x1 = r;end16.mulFastDown用最速下降法求非线性方程组的一组解function [r,m]=mulFastDown(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endlamda = fx/sum(diag(transpose(J)*J));r=x0-J*lamda; %核心迭代公式fr = subs(F,findsym(F),r);tol=dot(fr,fr);x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;17.mulGSND用高斯牛顿法求非线性方程组的一组解function [r,m]=mulGSND(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endDF = inv(transpose(J)*J)*transpose(J);r=x0-DF*fx; %核心迭代公式tol=norm(r-x0);x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;18.mulConj用共轭梯度法求非线性方程组的一组解function [r,m]=mulConj(F,x0,h,eps)format long;if nargin==3eps=1.0e-6;endn = length(x0);x0 = transpose(x0);fx0 = subs(F,findsym(F),x0);p0 = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)*(1+h);p0(:,i) = -(subs(F,findsym(F),x1)-fx0)/h;endm=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endlamda = fx/sum(diag(transpose(J)*J));r=x0+p0*lamda; %核心迭代公式fr = subs(F,findsym(F),r);Jnext = zeros(n,n);for i=1:nx1 = r;x1(i) = x1(i)+h;Jnext(:,i) = (subs(F,findsym(F),x1)-fr)/h;endabs1 = transpose(Jnext)*Jnext;abs2 = transpose(J)*J;v = abs1/abs2;if (abs(det(v)) < 1)p1 = -Jnext+p0*v;elsep1 = -Jnext;endtol=norm(r-x0);p0 = p1;x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;19.mulDamp用阻尼最小二乘法求非线性方程组的一组解function [r,m]=mulDamp(F,x0,h,u,v,eps)format long;if nargin==5eps=1.0e-6;endFI = transpose(F)*F/2;n = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsj = 0;fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;afx = subs(F,findsym(F),x1);J(:,i) = (afx-fx)/h;endFIx = subs(FI,findsym(FI),x0);for i=1:nx2 = x0;x2(i) = x2(i)+h;gradFI(i,1) = (subs(FI,findsym(FI),x2)-FIx)/h;ends=0;while s==0A = transpose(J)*J+u*eye(n,n);p = -A\gradFI;r = x0 + p;FIr = subs(FI,findsym(FI),r);if FIr<FIxif j == 0u = u/v;j = 1;elses=1;endelseu = u*v;j = 1;if norm(r-x0)<epss=1;endendendx0 = r;tol = norm(p);m=m+1;if(m>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endendformat short;。
数学实验报告非线性方程求解一、实验目的1.掌握用 MATLAB 软件求解非线性方程和方程组的基本用法,并对结果作初步分析;2.练习用非线性方程和方程组建立实际问题的模型并进行求解。
二、实验内容题目1【问题描述】(Q1)小张夫妇以按揭方式贷款买了1套价值20万元的房子,首付了5万元,每月还款1000元,15年还清。
问贷款利率是多少?(Q2)某人欲贷款50 万元购房,他咨询了两家银行,第一家银行开出的条件是每月还4500元,15 年还清;第二家银行开出的条件是每年还45000 元,20 年还清。
从利率方面看,哪家银行较优惠(简单假设:年利率=月利率×12)?【分析与解】假设初始贷款金额为x0,贷款利率为p,每月还款金额为x,第i个月还完当月贷款后所欠银行的金额为x i,(i=1,2,3,......,n)。
由题意可知:x1=x0(1+p)−xx2=x0(1+p)2−x(1+p)−xx3=x0(1+p)3−x(1+p)2−x(1+p)−x……x n=x0(1+p)n−x(1+p)n−1−⋯−x(1+p)−x=x0(1+p)n−x (1+p)n−1p=0因而有:x0(1+p)n=x (1+p)n−1p (1)则可以根据上述方程描述的函数关系求解相应的变量。
(Q1)根据公式(1),可以得到以下方程:150p(1+p)180−(1+p)180+1=0设 f(p)=150p(1+p)180−(1+p)180+1,通过计算机程序绘制f(p)的图像以判断解p的大致区间,在Matlab中编程如下:for i = 1:25t = 0.0001*i;p(i) = t;f(i) = 150*t*(1+t).^180-(1+t).^180+1;end;plot(p,f),hold on,grid on;运行以上代码得到如下图像:f(p)~p关系曲线图通过观察上图可知p∈[0.002,0.0022]。
Solution1:对于p∈[0.002,0.0022],采用二分法求解,在Matlab 中编程如下:clear;clc;x0=150000;n=180;x=1000;p0=0.002;p1=0.0022;while (abs(p1-p0)>1e-8)f0=x0*(1+p0).^n+x*(1-(1+p0).^n)/p0;f1=x0*(1+p1).^n+x*(1-(1+p1).^n)/p1;p2=(p0+p1)/2;f2=x0*(1+p2).^n+x*(1-(1+p2).^n)/p2;if (f0*f2>0 && f1*f2<0)p0=p2;elsep1=p2;end;end;p0结果得到p0=0.00208116455078125=0.2081%.所以贷款利率是0.2081%。
在MATLAB中解线性方程组或非线性方程(如一阶常微分方程或某些类型的一阶偏微分方程)通常可以使用内建的函数或库来完成。
这里是一些例子。
一、解线性方程组:
假设你有一个方程组Ax=b,其中A是系数矩阵,b是右侧的常数向量。
在MATLAB中,你可以使用`inv()`和`solve()`函数来求解这个方程组。
```matlab
% 定义系数矩阵A和右侧常数向量b
A = [1, 2; 3, 4];
b = [5; 6];
% 使用inv()和solve()函数求解
x = solve(A, b);
```
二、解非线性方程:
MATLAB也提供了内置函数来求解非线性方程。
例如,你可以使用`fzero()`函数来求解方程f(x)=0。
这个函数会使用一种名为"Secant"的迭代方法来寻找根。
```matlab
% 定义一个非线性函数,例如f(x) = x^2 - 4
f = @(x) x^2 - 4;
% 使用fzero()函数求解方程
x = fzero(f, [1, 1]); % 这将寻找一个可能的解,但是不一定是最优解
```
请注意,MATLAB中这些函数可能不适用于所有类型的方程。
在某些情况下,你可能需要使用其他方法或库,如数值积分或数值微分等。
如果你需要解决特定类型的方程,请提供更多详细信息,以便我可以提供更具体的帮助。
一实验目的1. 掌握雅可比方法和高斯-赛德尔方法的基本思想。
2. 编程实现雅可比方法和高斯-赛德尔方法求解非线性方程组的根。
3. 进一步熟悉matlab软件的使用。
二实验内容1、用雅可比方法求方程组的根,要求精度为1e-4。
(matlab)l1 计算公式l2 算法分析把方程组写成矩阵相成的形式,Ax=b,经变换得到x(k)的值,按上面的公式得到x(k+1),再用它代替x(k)继续以上过程,直到精度达到要求为止,即。
l3 源程序function [xx,k]=fun3(err)err=1e-4;a=[8,-3,2;4,11,-1;6,3,12];b=[20;33;36];xx=zeros(3,1);m=zeros(3,1);L=zeros(3,1);D=zeros(3,1);k=0;while k<100xx=L;for i=1:3for j=1:3if i~=jm(j)=a(i,j)*xx(j);endendD=sum(m);m=0;L(i)=(b(i)-D)/a(i,i);endif max(abs(xx-L))<errbreak;endk=k+1;end[xx;k]2、用高斯-赛德尔方法计算的根,要求精确到小数点后四位。
l6 计算公式l7 算法分析适当提供迭代初值,按G-S迭代分量公式讲旧值加工成新值,若迭代偏差,则输出结果,否则执行下一步。
若迭代次数k尚未达到最大迭代次数,则继续迭代,否则输出失败标志,终止计算。
l8 源程序clear;a=[8,-3,2;4,11,-1;6,3,12];b=[20;33;36];[m,n]=size(a);err=1e-4;x=zeros(n,1);s=zeros(n,1);k=0;while k<100xk=x;for i=1:nfor j=1:nif i~=jL(j)=a(i,j)*x(j);endends=sum(L);L=0;x(i)=(b(i)-s)/a(i,i);endif max(abs(xk-x))<errbreak;endk=k+1;end[x;k]三本次实验总结雅可比方法的迭代程序简单,每迭代一次只需计算一次矩阵和向量的乘法,存在收敛性,高斯-赛德尔方法算出一个新的分量就把前一次的分量冲掉,只需一组工作单元,算法更简单,在雅可比方法收敛的情况下,高斯-赛德尔方法收敛速度比雅可比方法的收敛速度快。
数学实验报告
Matlab的简单应用
——非线性方程(组)求解
姓名班级学号学院
2013年5月12日
一、实验目的
1.熟悉MATLAB软件中非线性方程(组)的求解命令及其用法。
2.掌握求非线性方程近似根的常用数值方法——迭代法。
3.了解分叉与混沌概念。
二、实验问题
1.利用弦截法编程对方程x^5+x-1=0进行求解实验,并与二分法、牛顿切线法进行
比较;
2.方程f(x)=x^2+x-4=0在(0,4)内有唯一的实根,现构造以下三种迭代函数:
(1)g1(x)=4-x^2,迭代初值x0=4;
(2)g2(x)=4/(1+x),迭代初值x0=4;
(3)g3(x)=x-(x^2+x-4)/(2x+1),迭代初值x0=4;
分别用给出的3种迭代函数构造迭代数列x(k+1)=g1(x(k)),i=1,2,3,观察这些迭代数列是否收敛,若收敛能否收敛到方程f(x)=0的解。
除此之外,你还能构造出其他收敛的迭代吗?
4.分别取不同的参数值r,做迭代数列x(n+1)=rx(n)(1-x(n)),n=0,1,2……,观察分
叉与混沌现象。
步骤1:首先,分别取参数r为0,0.3,0.6,0.9,1.2,1.5,1.8,2.1,2.4,2.7, 3.0,3.3,
3.6,3.9等14个值,按迭代序列迭代150步,分别产生14个迭代序列
{x(k)},k=0,1,…,150;其次,分别取这14个迭代序列的后50个迭代值
(x100,x101,…,x150),画在以r为横坐标的同一坐标面rox上,每一个r取值对应的迭代值点为一列。
步骤2:对(1)中图进行观察分析,容易发现:
(1)当r为0,0.3,0.6,0.9,1.2,1.5,1.8,2.1,2.4,2.7时,每个r对应的50个迭代值凝聚在一点,这说明对这些r的取值所产生的迭代序列是收敛的。
(2)当r为3,3.3时,r对应的50个迭代值凝聚在两个点,这说明这些r值所对应的迭代序列不收敛,但凝聚在两个点附近;同时也说明当r在2.7和3之间取值时,对
应的迭代序列从收敛到不收敛,轨道由一只分为两支开始出现分叉现象。
(3)当r由3.3到3.6再到3.9越来越大时,对应的50个迭代值凝聚的点也越来越多,表明r对应的迭代序列变化情况逐渐复杂,轨道分岔也越来越多,但会不
会还是按照一支分叉为两支的变化规律来变化呢?
步骤3:为了进一步研究上面所提到的问题,现在对r在2.7到3.9之间的取值进行加密迭代并作图,取步长为0.005时得到图像。
实验过程及结果分析
1.代码如下:
f=inline('x^5+x-1');
II/ 11
a=0;b=1;
x0=a;
k=1;
while abs(f(x0))>1.0e-5
x1=b-(b-x0)*f(b)/(f(b)-f(x0));
x0=x1;
vpa(x0,7)
k=k+1
plot(k,x0,'*')
hold on
end
结果:
……
k =
9
ans =
0.7548042
k =
10
ans =
0.7548515
k =
11
ans =
0.7548683
k =
12
ans =
0.7548743
k =
13
III
附:牛顿切线法
f=inline('x^5+x-1');
d1f=inline('5*x^4+1');
d2f=inline('20*x^3');
a=0;b=1;
if f(a)*d2f(a)>0
x0=a;
else
x0=b;
end
k=1;
while abs(f(x0))>0.00005
x1=x0-f(x0)/d1f(x0);
k=k+1
x0=x1
end
结果:
k =
2
x0 =
0.8333
IV/ 11
k =
3
x0 =
0.7644
k =
4
x0 =
0.7550
k =
5
x0 =
0.7549
2. (1)
x=0:0.01:4;
y1=x;
y2=4./(1+x);
plot(x,y1,x,y2)
hold on
x0=4;x1=4./(1+x0);s=[];ss=[];k=1;
while abs(x1-x0)>1.0e-5
x0=x1;s=[s,x0];
x1=4/(1+x0);ss=[ss,x1];
plot(s,ss,'r.')
pause(0.5)
vpa(x0,7)
k=k+1
end
结果:
……
k =
24
ans =
1.561563
k =
25
ans =
1.561546
k =
26
ans =
1.561557
V
k =
27
2.(2)
x=0:0.1:4;
y1=x;
y2=x-(x.^2+x-4)./(2.*x+1);
plot(x,y1,x,y2)
hold on
x0=4;x1=x0-(x0.^2+x0-4)./(2.*x0+1);s=[];ss=[];
while abs(x1-x0)>1.0e-5
x0=x1;s=[s,x0];
x1=x0-(x0.^2+x0-4)./(2.*x0+1);ss=[ss,x1];
plot(s,ss,'r*')
pause(0.5)
k=k+1
vpa(x1,7)
end
结果:
ans =
2.222222
k =
VI/ 11
2
ans =
1.641723
k =
3
ans =
1.563053
k =
4
ans =
1.561553
k =
5
2.(3)
x=0:0.01:4;
y1=x;
y2=4./(1+x);
plot(x,y1,x,y2)
hold on
f=inline('x^2+x-4');
d1f=inline('2*x+1');
x0=4;x1=4./(1+x0);s=[];ss=[];
VII
while abs(x1-x0)>1.0e-5
x0=x1;s=[s,x0];
x1=x0-f(x0)/d1f(x0);ss=[ss,x1];
plot(s,ss,'r.')
pause(0.5)
vpa(x1,7)
k=k+1
end
结果:
ans =
0.8
k =
2
ans =
1.784615
k =
3
ans =
1.572442
k =
4
ans =
1.561581
k =
5
ans =
1.561553
k =
6
VIII/ 11
3.
for r=0:0.3:3.9
n=0;
x0=0.1;
while n<=150
x1=r*x0*(1-x0);
x0=x1;
n=n+1
if n>50
plot(r,x0,'*')
hold on
end
end
hold on
end
结果:
IX
3'.
for r=2.7:0.005:3.9
n=0;
x0=0.1;
while n<=150
x1=r*x0*(1-x0);
x0=x1;
n=n+1
if n>50
plot(r,x0,'.')
hold on
end
end
hold on
end
结果:
X/ 11
西安交通大学数学实验报告
三、实验总结与体会
1.通过本次数学实验,我们基本掌握了利用matlab软件求解非线性方程组的
方法,这在我们日常学习中有很重要的应用,将我们从计算中解放出来从而专注于数学思想与方法的研究。
2.基本学会了利用弦截法、二分法、牛顿切线法求解非线性方程组的解并对这
三种方法进行了比较,三种方法分别适用于不同实际情况。
3.利用不同迭代法求解求非线性方程组的近似解,其中曾因为未考虑收敛关系
导致无法求解,最后发现问题并选取合适的收敛数列得出非线性方程组的近似解。
4.我们通过取加密r值不同的值得到对应的迭代序列,直观地看出迭代序列的
变化趋势越来越复杂,形象生动地认识了分叉与混沌现象。
XI。