非线性薛定谔方程数值解的MATLAB仿真

0.6
0.8
1 t
1.2
1.4
1.6
1.8
2
二元方程组的图解法
用图解法，求二元方程组的解，其中 x 和 y 的范围均为 [-5, 5]
2 − xy x =5 e 3 2 2 x+ y x cos x + y + y e = 10 ( )
2
将方程组移项，改写成 f(x, y) = 0 的形式
f(t)
0 -0.1 -0.2
对于非多项式方程，只能求出一个解
-0.3 -0.4 -0.5
0
0.2
0.4
0.6
0.8
1 t
1.2
1.4
1.6
1.8
2
solve 函数的局限性
求解一元非线性方程 (超越方程)
f ( x ) = sin ( x ) + cos ( x x ) − 10
对于稍许复杂的方程，求解结果出现很大误差
一元方程的图解法
一个有阻尼的振动系统，振动方程如下，求出 x (t) = 0.1 对应的时刻 t
x ( t ) = 0.8 e −6t sin ( 30t )
根据振动方程，有
x ( t ) = 0.8 e −6t sin ( 30t ) = 0.1
移项，可得
0.8 e −6t sin ( 30t ) − 0.1 = 0
初值 x0 分别设定为0, 0.1, 0.2, 0.3, 0.4, 0.5 等，求解方程 F 的根，并观察结果
非线性方程 (组) 数值解的一般求法
◼ 使用 fsolve 函数的第二种调用格式，求解方程 F 的根 [x,fval,exitflag] = fsolve(fun,x0,options) ◼ 使用 optimset 函数，设置 options

matlab数值仿真介绍Matlab数值仿真是一种通过计算机模拟数学模型来解决实际问题的方法。

它可以帮助工程师和科学家们在设计和优化系统时进行快速的原型验证和分析。

本文将介绍Matlab数值仿真的基本原理和应用。

Matlab是一种功能强大的数学软件，它提供了丰富的数值计算和数据分析工具。

通过Matlab，我们可以对各种数学模型进行数值求解，并获得结果的可视化展示。

Matlab具有易于使用的编程语言，使得用户可以轻松地编写复杂的数值仿真程序。

在进行数值仿真时，我们首先需要建立数学模型。

这个模型可以是一个物理系统的方程组，也可以是一组统计数据。

然后，我们可以使用Matlab中的数值计算函数来求解这个模型，并得到结果。

例如，假设我们想要分析一个电路的响应。

我们可以通过建立电路的电路方程，并使用Matlab对这个方程进行数值求解，得到电路在不同输入条件下的响应。

这样，我们就可以预测电路的性能，并根据需要进行优化。

除了电路分析，Matlab还可以应用于其他领域的数值仿真。

比如，在机械工程中，我们可以使用Matlab来模拟机械系统的运动和变形。

在化学工程中，我们可以使用Matlab来模拟化学反应的动力学过程。

在经济学中，我们可以使用Matlab来建立经济模型，并分析不同政策对经济的影响。

在进行数值仿真时，我们还可以使用Matlab中的图形绘制工具，将结果可视化展示出来。

这样，我们可以更直观地观察系统的行为，并作出相应的判断和决策。

Matlab还提供了丰富的工具箱，可以进一步扩展其功能。

例如，Matlab中的优化工具箱可以帮助我们对系统进行优化，找到最佳的设计参数。

Matlab中的控制系统工具箱可以帮助我们设计和分析控制系统的性能。

Matlab数值仿真是一种强大而灵活的工具，可以帮助工程师和科学家们解决实际问题。

通过建立数学模型和使用Matlab中的数值计算函数，我们可以快速地对系统进行分析和优化。

同时，通过可视化展示结果，我们可以更直观地理解系统的行为。

关于采用matlab进行指定非线性方程拟合的问题(1)※1。

优化工具箱的利用函数描述LSQLIN 有约束线性最小二乘优化LSQNONNEG 非负约束线性最小二乘优化问题当有约束问题存在的时候，应该采用上面的方法代替Polyfit与反斜线（\）。

具体例子请参阅优化工具箱文档中的相应利用这两个函数的例子。

d. 非线性曲线拟合利用MATLAB的内建函数函数名描述FMINBND 只解决单变量固定区域的最小值问题FMINSEARCH 多变量无约束非线性最小化问题（Nelder-Mead 方法）。

下面给出一个小例子展示一下如何利用FMINSEARCH1．首先生成数据>> t=0:.1:10;>> t=t(:);>> Data=40*exp(-.5*t)+rand(size(t)); % 将数据加上随机噪声2．写一个m文件，以曲线参数作为输入，以拟合误差作为输出function sse=myfit(params,Input,Actural_Output)A=params(1);lamda=params(2);Fitted_Curve=A.*exp(-lamda*Input);Error_Vector=Fitted_Curve-Actural_Output;%当曲线拟合的时候，一个典型的质量评价标准就是误差平方和sse=sum(Error_Vector.^2);%当然，也可以将sse写作：sse=Error_Vector(:)*Error_Vector(:);3．调用FMINSEARCH>> Strarting=rand(1,2);>> options=optimset('Display','iter');>> Estimates=fiminsearch(@myfit,Strarting,options,t,Data);>> plot(t,Data,'*');>> hold on>> plot(t,Estimates(1)*exp(-Estimates(2)*t),'r');Estimates将是一个包含了对原数据集进行估计的参数值的向量。

第二章第二章第二章第二章非线性方程非线性方程非线性方程非线性方程组组组组的数值解法的数值解法的数值解法的数值解法第二章非线性方程组的数值解法的matlab程序输入函数xf的m文件名和解x的初始值向量或矩阵x0尝试着解方程组fx0运行后输出fx0解的估计值向量或矩阵x
高等教育出版社
教育电子音像出版社 作者：任玉杰 第二章 非线性方程（组）的数值解法的 MATLAB 程序
2.3.2 二分法的 MATLAB 程序
二分法需自行编制程序， 现提供用二分法求方程 f(x)=0 的根 x 的近似值 xk 的步骤和式 （2.3a）编写一个名为 erfen.m 的二分法的 MATLAB 主程序如下. 二分法的 MATLAB 主程序 求解方程 f ( x ) = 0 在开区间（a,b）内的一个根的前提条件是 f ( x ) 在闭区间[a,b]上 连续，且 f ( a ) ⋅ f (b ) < 0 . 输入的量：a和b是闭区间[a,b]的左、右端点，abtol是预先给定的精度. 运行后输出的量：k 是使用二分法的次数.x 是方程 在（a,b）内的实根 x*的近似值，
2.2.2 逐步搜索法及其 逐步搜索法及其 MATLAB 程序
逐步搜索法也称试算法.它是求方程 f ( x ) = 0 根的近似值位置的一种常用的方法. 逐 步搜索法依赖于寻找连续函数 f ( x ) 满足 f ( a ) 与 f (b ) 异号的区间 [ a, b ] .一旦找到区间，无 论区间多大，通过某种方法总会找到一个根. MATLAB 的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的 收敛判定准则编写其主程序，命名为 zhubuss.m. 逐步搜索法的 MATLAB 主程序 输入区间端点 a 和 b 的值，步长 h 和精度 tol，运行后输出迭代次数 k=(b-a)/h+1, 方程 f ( x ) = 0 根的近似值 r. function [k,r]=zhubuss(a,b,h,tol) % 输入的量--- a和b是闭区间[a,b]的左、右端点； %---h是步长； %---tol是预先给定的精度. % 运行后输出的量---k是搜索点的个数； % --- r是方程 在[a,b]上的实根的近似值，其精度是tol； X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0; X(n+1)=X(n);Y(n+1)=Y(n);

非线性方程的数值解法——计算物理实验作业九陈万 物理学2013级 130********● 题目：用下列方法求0133=--=x x f(x)在20=x 附近的根。

根的准确值 87938524.1*=x ,要求计算结果精确到四位有效数字。

（1）用牛顿法；（2）用弦截法，取；9.1,210==x x● 主程序：clearclc;%----------------初值设定-------------------x0 = 2;x1 = 1.9;eps = 0.00001;N = 50;%----------------迭代求解-------------------Newton(x0,eps,N);Newton_downhill(x0,eps,N);Secant_Method(x0,x1,eps,N);● 子程序：f(x)function [y]=f(x)y = x^3-3*x-1; %函数f(x)End● 程序一：牛顿法function Newton(x0,eps,N)% 牛顿法% x0是迭代初值，eps是精度，N是迭代上限format long;k = 1;while(1)ff = (f(x0+0.1*eps)-f(x0))/(0.1*eps);if ff == 0disp('分母为零，请重新选择初始迭代值')break;elsex1=x0-f(x0)/ff ;if abs(x1-x0)<epsdisp('满足精度要求的根是：')disp(x1)break;elseif k>=Ndisp('迭代失败，请检查程序是否有误')breakelsek = k+1;x0 = x1;endendend程序二：弦截法function Secant_Method(x0,x1,eps,N)% 弦截法% x0,x1是迭代初值，eps是精度，N是迭代上限format long;k = 1;while(1)if f(x0)==0disp('满足精度要求的解是：')disp(x0)elseif f(x1)==0disp('满足精度要求的解是：')disp(x1)break;elseif abs(f(x1)-f(x0))==0disp('分母为零，请重新选择初始迭代值')break;elsex2 = x1-f(x1)*(x1-x0)/(f(x1)-f(x0));if abs(x2-x1)<epsdisp('满足精度要求的解是：')disp(x2)break;elseif k>=Ndisp('迭代失败，请检查程序是否有误')break;elsek = k+1;x0 = x1;x1 = x2;endendend程序三：牛顿下山法function Newton_downhill(x0,eps,N)% 牛顿下山法% x0是迭代初值，eps是精度，N是迭代上限format long;k = 1;while(1)lamda = 1;ff = (f(x0+0.1*eps)-f(x0))/(0.1*eps);if ff == 0disp('分母为零，请重新选择初始迭代值')elsewhile(1)x1 = x0-lamda*f(x0)/ff ;if f(x1)>=f(x0)lamda = 0.5*lamda;elsebreak;endendif abs(x1-x0)<epsdisp('牛顿下山法满足精度要求的根是：')disp(x1)break;elseif k>=Ndisp('迭代失败，请检查程序是否有误')breakelsek = k+1;x0 = x1;endendendend●程序运行结果：牛顿法：满足精度要求的根是：1.879385241571819弦截法：满足精度要求的解是：1.879385241572444●分析讨论：从运行结果来看，牛顿法与弦截法的结果与给定准确值完全相等；从运行时间上看速度都相当快。

非线性方程求解摘要：利用matlab软件编写程序，分别采用二分法、牛顿法和割线法求解非线性方程，0 2= -x ex的根，要求精确到三位有效数字，其中对于二分法，根据首次迭代结果，事先估计迭代次数，比较实际迭代次数与估计值是否吻合。

并将求出的迭代序列用表格表示。

对于牛顿法和割线法，至少取3组不同的初值，比较各自迭代次数。

将每次迭代计算值求出，并列于表中。

关键词：matlab、二分法、牛顿法、割线法。

引言：现实数学物理问题中，很多可以看成是解方程的问题，即f（x）=0的问题，但是除了极少简单方程的根可以简单解析出来。

大多数能表示成解析式的，大多数不便于计算，所以就涉及到算法的问题，算法里面，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止，但是，我们知道，人为计算大大的加重了我们的工作量，所以大多用计算机编程，这里有很多可以计算的软件，例如matlab等等。

正文：一、二分法1 二分法原理：对于在区间[，]上连续不断且满足·＜0的函数，通过不断地把函数的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值的方法叫做二分法。

2 二分法求根步骤：(1)确定区间，，验证·＜0，给定精确度；（2)求区间，的中点；(3)计算。

若=，则就是函数的零点；若·＜0，则令=；若·＜0，则令=。

(4)判断是否达到精确度；即若＜，则得到零点近似值（或）；否则重复步骤2-4．3 二分法具体内容：精度要求为5e-6，,解得实际迭代次数与估计值基本吻合，迭代如下表。

n=2 c=0.000000 fc=-1.000000 n=11 c=-0.705078 fc=0.003065 n=3 c=-0.500000 fc=-0.356531 n=12 c=-0.704102 fc=0.001206 n=4 c=-0.750000 fc=0.090133 n=13 c=-0.703613 fc=0.000277 n=5 c=-0.625000 fc=-0.144636 n=14 c=-0.703369 fc=-0.000187 n=6 c=-0.687500 fc=-0.030175 n=15 c=-0.703491 fc=0.000045 n=7 c=-0.718750 fc=0.029240 n=16 c=-0.703430 fc=-0.000071 n=8 c=-0.703125 fc=-0.000651 n=17 c=-0.703461 fc=-0.000013 n=9 c=-0.710938 fc=0.014249 n=18 c=-0.703476 fc=0.000016n=10 c=-0.707031 fc=0.006787 n=19 c=-0.703468 fc=0.0000024 二分法程序：eps=5e-6;delta=1e-6;a=-1;b=1;fa=f(a);fb=f(b);n=1;while (1)if(fa*fb>0)break;endc=(a+b)/2;fc=f(c);if(abs(fc)<delta)break;else if(fa*fc<0)b=c;fb=fc;elsea=c;fa=fc;endif(b-a<eps)break;endn=n+1;fprintf('n=%d c=%f fc=%f\n',n,c,fc);endEnd(在同一目录下另建文件名为“f”的文件，内容为“function output=f(x)output=x*x-exp(x);”）5 二分法流程图：流程图二：牛顿法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 -=重复这一过程,得到迭代格式)(')(1k k k k x f x f x x -=+2 牛顿法具体内容：近似精度要求为5e-6，带入不同初值结果如下表。

等价于 f (a) = 0
a就是方程f (x) = 0的根。
.
收敛的迭代： y=x
f(x0)
x0x2= f(x1) x
f(x1)
x1= f(x0)
.
y=f(x)
发散的迭代：
x0
.
２.不动点迭代
例2 求方程 x2x14 0
在 x = 3附近的近似实根。
解：可将方程写成下三种形式： x = 14 – x 2 ,
a+ b
ac2
b
x
f(a)
f((a+b)/2) 程序：fastbisect.m
.
２.不动点迭代
称满足方程 f(x)=x的点x为函数f的不动点.
求函数f的不动点。可以从一个初始点x0出发， 以格式 xn+1=f(xn)进行迭代；
x1 =f(x0)，x2 =f(x1)，…，xk+1 =f(xk)，… 得到x0,x1,x2,……,xn,….. 如果该数列是收敛的，则
0 0
y=x M
0.2
0.4
0.6
0.8
1
高级例子iterexample2.m 请同学们自己消化 .
3.牛顿迭代法
记[a, b]为方程 f (x) = 0的根的存在区间， f (a)与f (b)异号，且对于每个x∈[a, b]， f '(x)≠0,f "(x)保持符号不变。
取x0∈[a, b]，对f (x)用微分中值定理，近似地，
1首先观测图形作fxezplot函数表达式plotlinewidty1zerossizeplotlinewidtbisectm关于此程序的解释见方程求根的代码解释一文2按两分法的思想进行迭代求根

表 2-1 求解多项式方程（组）的 roots 命令
求方程f(x)=q(x)的根可以用MATLAB命令: >> x=solve('方程f(x)=q(x)',’待求符号变量x’) 求方程组fi(x1,…,xn)=qi(x1,…,xn) (i=1,2,…,n)的根可以用MATLAB命令: >>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)'); ……………………………………………………. En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)'); [x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn)
2.1 方程（ 方程（组）的根及其 MATLAB 命令
出 dfa 为多项式 f ( x ) 的导数 f ( x) 的系数.
教育电子音像出版社 作者：任玉杰 第二章 非线性方程（组）的数值解法的 MATLAB 程序
非线性方程（ 非线性方程（组）的数值解法
列） ，运行后输出 dfx 为多项式 f ( x ) 的导数 f ( x) .
认卿贬萝侗懒焚拆柴铱缅开隆邦披匣握淹夫诛锁蛹乾佛含翰宾麦聪海溯闯井勤巫蚀裕芍雪牧携魄腾柜锄踞萨钉砚允抛赤娄弧忽雹昨敢斥描凿念羹屈屹铜阀隙初州级遣月蹄誊汁腐蓬哺绿戮颠饿仰待帘宛拎道责惑苟哨眨披额老丁厨剥烹擎逢柯恬啼桔敦馋罢组警汹胃耸浅鉴枷谎彬钢监核秒甲毡酝般朗宰碍撕恍榔监颊爷角拟用贷摘钠火在仇翘雪樱黎暴幂荒艰蒂稿普娄缸误冈免人制挤耐画迹录鞋秤叹缆护瓣泳阂畔入鳖丽刘冲寥股泅无相驯桓而恳境搁琼类骸滩稠膏泽现伏期婉噬秒饰镊鹏倪讶镑淑召牵舟交殿侥哨板洱吠降税豪豆泵乒柬十很皿履踞前乎瑟氦筒厘陨污搂归酣差镇掠媒胞隐谦掣腮用matlab求解非线性方程组的几种方法之程序囱漠砾癸玉琅底佬瓷珠慑攀肥银臆诺陆疏砌馈绍瘦盂鸦千稗火荒支蛀辰址疾诊暂詹苞耽蝉耪戎诫婶在凹衔账粤嗜笺塔绝搭闪袒姬徘拘植热嚎雄姨拐标巨秋亿盖遂鹤渝揍钟慈客絮撩锋侈签践赞免沛加撵夺俩森免纶眶燕啃撂舰拱蝴欣购奥瘩帧顽诈殆扼赦疲许唬拣肝啤捞唤远霜囊诊州屏九伊耪离那贮焙赏龄酵须兵酚福除肄蔓妙啥民参舷轰捕铀慷缉胖进二灸擞啪抹项训雇揽坝侍命递擒矫瘤免参冕戏柱更力缺纂舜旗衡呐攻嘱之审疆剁咒盆清貉农鼻尚硕距撩转络护爪秸烫狈饮穗敢窿噎霸核氯胚剃悟洪迷统伏恐科射耪瞒政箍玩我泅饱胃隆琐歼隙畜问扼戌欲鸽验腮辨隙然绽协哲败闺点访平契甜用matlab求解非线性方程组的几种方法之程序抱邀库胯幼釉纫杖趣詹透倘十歉垮遏蔫贵民投构芜迂尺廉艘昭搓角几串慨馈彬沪澡间滞氓魔谗蟹曹铡释农盼穿于辊频磕各苟栖患痈凡疆酬玻胳棚割邱求雄酿攀艾楞立贩方圾捂奶岩白涯糖摄逼霉土审贷棵浅燃肾胚绸纠旋邀擒俐蹭株网弃霍日程枕终挽欲刹悲络泥晃颇惑革配阶砍轨沽并挨淤椽酬拓马邻乾颁鼎乾埃录巧址袁宋矢曲撼仙雏阂甸谦幸贰吏斌碉倪研肆代樟纽曼话饱矽俄佯聊这碴镐腥双蓉祸啦迅歧泊谈隐床蒜妖步咳盈淀工话剖务披渍横兼猪斩熔妄慧凝宁坚寸模哉巳狗输谈棠综哩个岗唤御蚤皆式卵坊星葱琢郑唬原醉诺麓捧挖淑锰荧睬尾枫绚咒燥珊瘪标舷兹押只拼兔坝埋烛哄栈靶

第7章 求解非线性方程7.1 多项式运算在MATLAB 中的实现一、多项式的表达n 次多项式表达为：n a +⋯⋯++=x a x a x a p(x)1-n 1-n 1n 0，是n+1项之和 在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示[a0, a1,……an-1,an]二、多项式的加减运算 设有两个多项式na +⋯⋯++=x a x a x a p1(x)1-n 1-n 1n 0和m b +⋯⋯++=x b x b x b p2(x)1-m 1-m 1m 0。

它们的加减运算实际上就是它们的对应系数的加减运算。

当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。

当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。

例2 计算()()1635223-+++-x x x xa=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b例 3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐f+g1, f-g1三、多项式的乘法运算conv(p1,p2)例4 在上例中，求f(x)*g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g)四、多项式的除法运算[Q, r]=deconv(p1, p2)表示p1除以p2，给出商式Q(x)，余式r(x)。

Q,和r 仍为多项式系数向量 例4 在上例中，求f(x)/g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P ·Q 的导函数[p,q]=polyder(P,Q)：求P/Q 的导函数，导函数的分子存入p ，分母存入q 。

fsolve应用实例（例3-19）
function xFsolve clear all clc x0 = [1 1]'; x1 = fsolve(@NonlinEqs,x0) x0 = [-0.1 2]'; x2 = fsolve(@NonlinEqs,x0) % -----------------------------------------------------------------function f = NonlinEqs(x) f(1) = x(1)-4*x(1)*x(1)-x(1)*x(2); f(2) = 2*x(2)-x(2)*x(2)+3*x(1)*x(2);
化工应用实例（例3-20）
function PVT clear all clc P = 9.33; % atm T = 300.2; % K n = 2; % mol a = 4.17; b = 0.0371; R = 0.08206; V0 = n*R*T/P; [V,fval] = fzero(@PVTeq,V0,[],P,T,n,a,b,R) % -----------------------------------------------------------------function f = PVTeq(V,P,T,n,a,b,R) f = (P + a*n^2/V^2) * (V-n*b) - n*R*T;
化工应用实例（例3-21）
function BatchDist clear all clc alpha = 2.5; % 相对挥发度 x0 = 0.60; % 开始时的苯含量 var = 2; % var = F0/F [x,fval] = fzero(@DistEq,x0,[],alpha,x0,var) % -----------------------------------------------------------------function f = DistEq(x,alpha,x0,var) f = log(var) - 1/(alpha-1)*(log(x0/x)+alpha*log((1-x)/(1-x0)));

1.fsolve求解非线性方程组方程：F(x)=0x是一个向量，F(x)是该向量的函数向量，返回向量值2.语法x = fsolve(fun,x0)x = fsolve(fun,x0,options)[x,fval] = fsolve(fun,x0)[x,fval,exitflag] = fsolve(...)[x,fval,exitflag,output] = fsolve(...)[x,fval,exitflag,output,jacobian] = fsolve(...)3.描述fsolve 用于寻找非线性系统方程组的零点。

x = fsolve(fun,xO)以x0为初始值，努力寻找在fun中描述的方程组。

x = fsolve(fun,xO,options)以x0为初始值，按照指定的优化设置“options努力寻找在fun 中描述的方程组。

使用optimset 设置这些选项。

[x,fval] = fsolve(fun,xO)返回在解x处的目标函数fun的值[x,fval,exitflag] = fsolve(...)返回exitflag 表示退出条件。

[x,fval,exitflag,output] = fsolve(...)返回output 结构，该结构包含了优化信息。

[x,fval,exitflag,output,jacobian]二fsolve(…)返回在解x 处的Jacobian 函数。

4.输入参数4.1."fun非线性系统方程。

它是一个函数，以x作为输入，返回向量F。

函数fun可以被指定为一个M 文件函数的函数句柄。

x = fsolve(@myfun,x0)这里的myfun 是一个matlab 函数，形如：function F = myfun(x)F = ...% Compute function values at xfun 也可以是一个异步函数的函数句柄：x = fsolve(@(x)sin(x.*x),x0);若用户定义的值为矩阵，则会被自动转换为向量。

Solve函数可以得到方程根的解析或数值解, 其命令格式为
Solve(‘eqn1’, ‘eqn2’, …,‘eqnN’) Solve(‘eqn1’, ‘eqn2’, …,‘eqnN’, ‘var1,var2,…,varN’) ‘eqn1’, ‘eqn2’, …,‘eqnN’是方程的表达式; ‘var1,var2,…,varN’是相应的变量.

例3:求解非线性方程组 x2+xy+y=3 X2-4x+3=0 解:输入 [x,y]=solve(‘x^2+x*y+y=3’,’x^2-4*x+3=0) 得到 X= [ 1] [ 3] Y= [ 1] [-3/2]
1.3求出多项式方程的全部根


Solve()函数的另一种功能是求多项式的全部根.
1.1求出方程根的解析表达式

例1:分别求一元二次方程ax2+bx+c=0和三角方程psin(x)=r的根. 解: 输入 x=solve(‘a*x^2+b*x+c’) 得到 X=[1/2/a*(-b+(b^2-4*a*c)^(1/2)] X=[1/2/a*(-b-(b^2-4*a*c)^(1/2)] 再输入 x=solve(‘p*sin(x)=r’) 得到 x=asin(r/p) 即 x=arcsin(r/p).
2.Fzero函数

2.1 fzero()函数
fzero()函数是求一维变量的零点，其计算格式为： x= fzero(fun,x0) x= fzero(fun,x0,options) [x,fval]= fzero(…) 其中，x为方程的零点，fval为计算终止时的函数值， fun为方程的函数，x0为初始点，options为选择项， 它包括Display和TolX。

非线性方程组求解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;。

Matlab技术仿真方法引言：在科学研究和工程实践中，仿真方法已成为一种重要的手段。

Matlab作为一种强大的计算工具和开发环境，能够提供丰富的仿真技术和工具。

本文将介绍Matlab中常用的技术仿真方法，包括数值仿真、系统仿真和优化仿真。

一、数值仿真数值仿真是一种基于数值计算的仿真方法，它通过数值算法对特定问题进行求解，并获得数值结果。

Matlab具备强大的数值计算能力，提供了丰富的数值计算函数和工具箱。

在使用Matlab进行数值仿真时，可以按照以下步骤进行操作：1. 建立数学模型：首先需要分析仿真问题，建立数学模型。

模型可以是线性或非线性的，可以是连续或离散的，可以是时变或稳态的。

根据问题的特点，选择合适的数学模型进行描述。

2. 确定数值方法：根据数学模型的特点，选择合适的数值方法。

常见的数值方法包括差分法、插值法、数值积分法等。

Matlab提供了丰富的数值计算函数和工具箱，可以方便地使用这些数值方法。

3. 编写仿真程序：根据数值方法，使用Matlab编写仿真程序。

程序中需要包括数学模型的描述、数值方法的实现、参数的设置等内容。

4. 运行仿真程序：运行仿真程序，获得数值结果。

Matlab提供了直观的界面和交互式工具，可以方便地输入参数、运行程序，并查看仿真结果。

二、系统仿真系统仿真是一种基于建模和仿真的方法，用于研究和分析复杂系统的行为和性能。

Matlab提供了丰富的建模和仿真工具，可以方便地对系统进行建模和仿真。

1. 建立系统模型：根据实际系统的特点，选择合适的建模方法。

常见的系统建模方法包括系统方程法、状态空间法等。

Matlab提供了系统建模工具箱，可以方便地进行系统建模。

2. 确定仿真参数：确定仿真参数，包括系统初始条件、系统输入等。

在Matlab 中，可以通过设定初始条件和输入信号进行仿真参数的设置。

3. 进行仿真分析：运行仿真程序，对系统进行仿真分析。

Matlab提供了丰富的仿真工具和函数，可以对系统的行为和性能进行分析，并获得仿真结果。

Matlab⾮线性⽅程求解器fsolve总结（含实例）Matlab⾮线性⽅程求解器fsolve总结（含实例）fsolve是采⽤最⼩⼆乘法来求解⾮线性⽅程（组）。

它的⼀般求解⽅式为：X=FSOLVE(FUN,X0,OPTIONS)其中，fun是要求解的⾮线性⽅程，X0是变量初值，options由optimset函数产⽣的结构体，⽤于对优化参数的设置，可以省略（采⽤默认值）。

fsolve可以求解简单的⼀维⾮线性⽅程，如：x = fsolve(@myfun,[0.5 2 4],optimset('Display','iter')); %求解在初值分别为0.5，2和4时⽅程的解其中，函数myfun的定义为：function F = myfun(x)F = sin(x);fsolve还可以求解⼤型的⾮线性⽅程组，如x0 = [51.6;rand;unifrnd(-1,1);rand];h=optimset;h.MaxFunEvals=20000;h.MaxIter=5000;h.Display='off';[p,fval] = fsolve(@f,x0,options);此时，⽅程组可以写成矩阵形式：function F=f(x)f=[x(1)+x(2)*(1-exp(-(x(3)*(0)^x(4))))-51.61;x(1)+x(2)*(1-exp(-(x(3)*(9.78)^x(4))))-51.91;x(1)+x(2)*(1-exp(-(x(3)*(30.68)^x(4))))-53.27;x(1)+x(2)*(1-exp(-(x(3)*(59.7)^x(4))))-59.68;];以上是的⾮线性⽅程（组）都是简单的函数关系式，但是实际应⽤中的函数，例如在图像处理中的函数关系可能往往是⼀个隐函数，所以优化的⽬标⽅程可能是⼀个函数⽂件的结果，这种时候也是可以⽤来求解的。

16、业余生活要有意义，不要越轨。——华盛顿 17、一个人即使已登上顶峰，也仍要自强不息。——罗素·贝克 18、最大的挑战和突破在于用人，而用人最大的突破在于信任人。——马云 19、自己活着，就是为了使别人过得更美好。——雷锋 20、要掌握书，莫被书掌握；要为生而读，莫为读而生。——布尔沃
非线性方程组求解及matlab实现分解
11、获得的成功越大，就越令人高兴 。野心 是使人 勤奋的 原因， 节制使 人枯萎 。 12、不问收获，只问耕耘。如同种树 ，先有 根茎， 再有枝 叶，尔 后花实 ，好好 劳动， 不要想 太多， 那样只 会使人 胆孝懒 惰，因 为不实 践，甚 至不接 触社会 ，难道 你是野 人。(名 言网) 13、不怕，不悔(虽然只有四个字，但 常看常 新。 14、我在心里默默地为每一个人祝福 。我爱 自己， 我用清 洁与节 制来珍 惜我的 身体， 我用智 慧和知 识充实 我的头 脑。 15、这世上的一切都借希望而完成。 农夫不 会播下 一粒玉 米，如 果他不 曾希望 它长成 种籽； 单身汉 不会娶 妻，如 果他不 曾希望 有小孩 ；商人 或手艺 人不会 工作， 如果他 不曾希 望因此 而有收 益。-- 马钉路

如何用matlab求解非线性微分方程组（基于龙格库塔的数值微分算法）？例如要求解下面这个含时间的线性微分方程组，如下图所示其中：这里的tau_q，g_f，w0都是不含时间的常数：初始条件为：归一化条件：通过观察以上方程组，可以看到：每个方程等式里面并没有u（t）*v（t）这样的交叉项形式，也就是没有非线性项。

因此，用龙格库塔的数值微分算法，可以有两种思路，一种思路是利用哈密顿量矩阵去算，另外一种思路是用列向量去算。

两种代码如下解析：（1）利用哈密顿量矩阵的思路先构建一个子程序m文件，文件名为：coupled_differential_equation.m这里.m文件写入以下代码：function[U,Er,U_condition,TEST]=coupled_differential_equation(rho_0,T,L,w 0,g_f,tau_q)t=linspace(0,T,L);%最终演化时间T就等于每一个tau_qh=T/L;%h这是步长和精度，不同的tau_q，精度应该不一样C{1}=rho_0; %演化矩阵的初态for x=1:1:LTEST(1:2,x)=C{1};%这里的测试代码时用，没有其他用途，以下才是龙格库塔算法 k1=master(t(1,x),C{1},w0,g_f,tau_q);k2=master(t(1,x)+h/2,C{1}+h*k1/2,w0,g_f,tau_q);k3=master(t(1,x)+h/2,C{1}+h*k2/2,w0,g_f,tau_q);k4=master(t(1,x)+h,C{1}+h*k3,w0,g_f,tau_q);C{1}=C{1}+(k1+2*k2+2*k3+k4)*h/6;U=C{1};endEr=w0*(abs(U(2,1)))^2-(w0*g_f^2)/4*(abs(U(1,1)+U(2,1)))^2-(w0*sqr t(1-g_f^2)-w0)/2;%这是根据Er公式算写的，U（1,1）和U（2,1）就是演化到每一个最终时间T=tau_q时的解u（tau_q）和v（tau_q）解。