弦截法非线性方程求解2
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
fb=subs(sym(f),findsym(sym(f)),b);
root=a-(b-a)*fa/(fb-fa);%迭代初始值
while(tol>eps)
r1=root;
fx=subs(sym(f),findsym(sym(f)),r1);
s=fx*fa;
if(sΒιβλιοθήκη Baidu=0)
root=r1;
else
if(s>0)
%x(4)=dy/dt
mu=1/82.45;
lambda=1-mu;
r1=sqrt((x(1)+mu)^2+x(3)^2);
r2=sqrt((x(1)+lambda)^2+x(3)^2);
dx=[x(2);
2*x(4)+x(1)-lambda*(x(1)+mu)/r1^3-mu*(x(1)-lambda)/r2^3;
《MATLAB程序设计实践》课程考核
1、编程实现以下科学计算方法，并举一例应用之（参考书籍《精通MATLAB科学计算》，王正林著，电子工业出版社，2009年）“弦截法非线性方程求解”
算法说明：
（1）过两点 , 作一直线，它与 轴有一个交点，记为 ；
（2）如果f(a)f(x1)<0,过两点 ， 作一直线，它与 轴的交点记为 ，否则过两点 ， 作一直线，它与 轴的交点记为 ；
x(4);
-2*x(2)+x(3)-lambda*x(3)/r1^3-mu*x(3)/r2^3];
运行：apollorun.m
流程图：
结果：
但是此结果与实际的阿波罗卫星轨道不符，查找相关资料，发现题目出错
应改为
更正后结果：
与实际情况符合较好
关于为什么相对误差要取1e-6的原因：
如果相对误差设置为1e-3，即默认值，发现图形为：
>> r=Secant('sqrt(x)+log(x)-2',1,4)
输出计算结果为：
r =
1.8773
有计算结果可知， 在区间 上的一个根为： 。
2、分析单自由度阻尼系统的阻尼系数对其固有振动模态的影响。
算法分析：
单自由度阻尼系统振动方程：
其中k和c分别为弹簧的阻尼系数和黏性阻尼系数
弹簧阻力：
黏性阻力：
方法一：
算法说明：
根据该微分方程组可编写m文件，命名为apolloeq.m。在m文件编写中，由于出现了高阶微分方程，因此把这个高阶微分方程化为一阶微分方程组。然后调用ode45函数求解即可。
编写函数m文件：
functiondx=apolloeq(t,x)
%x(1)=x
%x(2)=dx/dt
%x(3)=y
root为求出的函数零点。
流程图：
源代码：
functionroot=Secant(f,a,b,eps)
%弦截法求函数在区间[a,b]上的一个零点
%函数名：f
%区间左端点：a
%区间右端点：b
%根的精度：eps
%求出函数的函数零点：root
if(nargin==3)
eps=1.0e-4;
end
f1=subs(sym(f),findsym(sym(f)),a);
模型：
结果：
（3）弱阻尼情况：即阻尼比 时
取c=4，k=0.5，m=2时：此时 2
模型：
结果：
三条图形的比较：
运行vibrationrun.m
得到结果：
其中：绿色为阻尼比等于一的情况，红色为阻尼比小于一的情况，蓝色为阻尼比大于一的情况。
流程图：
3、已知Apollo卫星的运动轨迹 满足下面方程：
其中， 试在初值 下求解，并绘制Apollo卫星轨迹图。
显然这条轨道不稳定，如果取1e-6就可以得到一条稳定的轨道，所以应该取1e-6
方法二：建立如下Simulink模型
其中子模块r1：
子模块r2：
子模块sub1：
子模块sub2：
子模块sub3：
子模块sub4：
设定相关参数：采用ODE45
当相对误差为1e-3时模拟得到图形：
发现精度太低 提高精度 把相对误差改成1e-6
root=b-(r1-b)*fb/(fx-fb);%用递推公式2
else
root=a-(r1-a)*fa/(fx-fa);%用递推公式1
end
end
tol=abs(root-r1);
end
end
例:弦截法求解非线性方程应用实例。采用弦截法求方程 在区间 上的一个根。
解：在MATLAB命令窗口中输入：
可解得：
根据微分方程建立Simulink模型：设仿真时间t=10s其中，利用To Workspace把模拟的数据输出到工作空间中，然后调用函数进行绘图和比较。
（1）欠阻尼情况：即阻尼比 时
取c=4，k=8，m=2时：此时 0.5
模型：
结果：
（2）临界阻尼情况：即阻尼比 时
取c=4，k=2，m=2时：此时 1
f2=subs(sym(f),findsym(sym(f)),b);
if(f1==0)
root=a;
end
if(f2==0)
root=b;
end
if(f1*f2>0)
disp('两端点函数值乘积大于0！');
return;
else
tol=1;
fa=subs(sym(f),findsym(sym(f)),a);
（3）如此下去，直到 就可以认为 为 在区间 上的一个根。
（4） 的递推公式为：
且 。
在MATLAB中编程实现的弦截法的函数为：Secant.
功能：用弦截法求函数在某个区间的一个零点。
调用格式：root=Secant(f,a,b,eps).
其中，f为函数名；
a为区间左端点；
b为区间左端点；
eps为根的精度；
模拟图形：
虽然此时曲线平滑，但是与实际的轨道不符，查阅相关资料发现原来题目中有条件出错。
应改为
此时，把子模块r2改为：
进行模拟后得到结果
此时与实际条件相符，此为正确结果。
root=a-(b-a)*fa/(fb-fa);%迭代初始值
while(tol>eps)
r1=root;
fx=subs(sym(f),findsym(sym(f)),r1);
s=fx*fa;
if(sΒιβλιοθήκη Baidu=0)
root=r1;
else
if(s>0)
%x(4)=dy/dt
mu=1/82.45;
lambda=1-mu;
r1=sqrt((x(1)+mu)^2+x(3)^2);
r2=sqrt((x(1)+lambda)^2+x(3)^2);
dx=[x(2);
2*x(4)+x(1)-lambda*(x(1)+mu)/r1^3-mu*(x(1)-lambda)/r2^3;
《MATLAB程序设计实践》课程考核
1、编程实现以下科学计算方法，并举一例应用之（参考书籍《精通MATLAB科学计算》，王正林著，电子工业出版社，2009年）“弦截法非线性方程求解”
算法说明：
（1）过两点 , 作一直线，它与 轴有一个交点，记为 ；
（2）如果f(a)f(x1)<0,过两点 ， 作一直线，它与 轴的交点记为 ，否则过两点 ， 作一直线，它与 轴的交点记为 ；
x(4);
-2*x(2)+x(3)-lambda*x(3)/r1^3-mu*x(3)/r2^3];
运行：apollorun.m
流程图：
结果：
但是此结果与实际的阿波罗卫星轨道不符，查找相关资料，发现题目出错
应改为
更正后结果：
与实际情况符合较好
关于为什么相对误差要取1e-6的原因：
如果相对误差设置为1e-3，即默认值，发现图形为：
>> r=Secant('sqrt(x)+log(x)-2',1,4)
输出计算结果为：
r =
1.8773
有计算结果可知， 在区间 上的一个根为： 。
2、分析单自由度阻尼系统的阻尼系数对其固有振动模态的影响。
算法分析：
单自由度阻尼系统振动方程：
其中k和c分别为弹簧的阻尼系数和黏性阻尼系数
弹簧阻力：
黏性阻力：
方法一：
算法说明：
根据该微分方程组可编写m文件，命名为apolloeq.m。在m文件编写中，由于出现了高阶微分方程，因此把这个高阶微分方程化为一阶微分方程组。然后调用ode45函数求解即可。
编写函数m文件：
functiondx=apolloeq(t,x)
%x(1)=x
%x(2)=dx/dt
%x(3)=y
root为求出的函数零点。
流程图：
源代码：
functionroot=Secant(f,a,b,eps)
%弦截法求函数在区间[a,b]上的一个零点
%函数名：f
%区间左端点：a
%区间右端点：b
%根的精度：eps
%求出函数的函数零点：root
if(nargin==3)
eps=1.0e-4;
end
f1=subs(sym(f),findsym(sym(f)),a);
模型：
结果：
（3）弱阻尼情况：即阻尼比 时
取c=4，k=0.5，m=2时：此时 2
模型：
结果：
三条图形的比较：
运行vibrationrun.m
得到结果：
其中：绿色为阻尼比等于一的情况，红色为阻尼比小于一的情况，蓝色为阻尼比大于一的情况。
流程图：
3、已知Apollo卫星的运动轨迹 满足下面方程：
其中， 试在初值 下求解，并绘制Apollo卫星轨迹图。
显然这条轨道不稳定，如果取1e-6就可以得到一条稳定的轨道，所以应该取1e-6
方法二：建立如下Simulink模型
其中子模块r1：
子模块r2：
子模块sub1：
子模块sub2：
子模块sub3：
子模块sub4：
设定相关参数：采用ODE45
当相对误差为1e-3时模拟得到图形：
发现精度太低 提高精度 把相对误差改成1e-6
root=b-(r1-b)*fb/(fx-fb);%用递推公式2
else
root=a-(r1-a)*fa/(fx-fa);%用递推公式1
end
end
tol=abs(root-r1);
end
end
例:弦截法求解非线性方程应用实例。采用弦截法求方程 在区间 上的一个根。
解：在MATLAB命令窗口中输入：
可解得：
根据微分方程建立Simulink模型：设仿真时间t=10s其中，利用To Workspace把模拟的数据输出到工作空间中，然后调用函数进行绘图和比较。
（1）欠阻尼情况：即阻尼比 时
取c=4，k=8，m=2时：此时 0.5
模型：
结果：
（2）临界阻尼情况：即阻尼比 时
取c=4，k=2，m=2时：此时 1
f2=subs(sym(f),findsym(sym(f)),b);
if(f1==0)
root=a;
end
if(f2==0)
root=b;
end
if(f1*f2>0)
disp('两端点函数值乘积大于0！');
return;
else
tol=1;
fa=subs(sym(f),findsym(sym(f)),a);
（3）如此下去，直到 就可以认为 为 在区间 上的一个根。
（4） 的递推公式为：
且 。
在MATLAB中编程实现的弦截法的函数为：Secant.
功能：用弦截法求函数在某个区间的一个零点。
调用格式：root=Secant(f,a,b,eps).
其中，f为函数名；
a为区间左端点；
b为区间左端点；
eps为根的精度；
模拟图形：
虽然此时曲线平滑，但是与实际的轨道不符，查阅相关资料发现原来题目中有条件出错。
应改为
此时，把子模块r2改为：
进行模拟后得到结果
此时与实际条件相符，此为正确结果。