非线性方程模型及求解实例

---

---

一、气体方程模型

例 1.
1摩尔理想气体的压强P，体积V，温度T满足关系PV=RT，其中常数R=0.08205（latm/Kmol）而对于实际气体这个关系需修正为
（P + $\frac{a}{V^c}$) (V - b) = RT
a,b,c为所给气体决定的常数。现已知a=12.87, b=0.1142, c=2，求气体在P=2(atm)，T=313(K)下的体积V。

分析
首先，根据题目中所给出的数值与公式方程可得出本题的基本数学模型，方程是一个非线性方程组，上篇博客写到，可以使用fzero函数解方程，当然，在解方程之前，可画出方程的大致图像，判断出根的大致范围
代码如下：

%气体的一些常数
a=12.87;b=0.1142;c=2;
R=0.08205;
%压强和温度
P=2;T=313;
%理想气体方程：P*V=R*T；
%实际气体方程：
%f(V) = (p+a / V ^ c) * (V - b) - R * T 
V0 = R * T / P;
f=@(V) (P + a ./  V .^ c) .* (V - b) - R * T;
fplot(f, [0.1, 15]);
V = fzero(f, V0)

可得到图像：

从图像中可以看出，方程的解有三个，但是我们只求近似解 V0 = R * T / P附近的根，结果为：

V =
   12.4426

二、房贷还款问题

例2. （购房贷款的利率）
下面是一则房产广告：

不难算出，你向银行总共借了25.2万，
30年内共要还51.696万，
这个案例中贷款年利率是多少呢？

分析
设$x_k$为第k个月末的欠款数， a为月还款数， r为月利率。
$x_k$ = (1+ r) $x_{k-1}$ - a
= (1+r)² $x_{k-2}$ – (1+r) a – a
=……
= (1+r)^k $x_0$ – a[1+(1+r)+……+(1+r)^k-1]
= (1+r)^k $x_0$ – a $\frac{(1+r{)}^k-1}{r}$
根据题意可知$x_{360}$ = 0，由此可以列出方程，画出图像求解即可

代码如下：

%房贷年利率问题
%设 xk 表示（从现在开始过了k个月后的）欠款（现在的时间为，k=0）
%设每个月还款a （单位：元），月利息为 r，则 x0 = 25.2
%x1 = x0（1 + r） - a （第一个月月初）
%x2 = x1（1 + r） - a （第二个月月初）
%……
%x360 = x359*（1 + r）- a（第360个月月初还 a）
%x360 = 0
x0=25.2;
a=0.1436;
f=@(r) x0*(1+r).^360 - a*( (1+r).^360 - 1 ) ./ r ;
r0=0.001;
r=fzero(f,r0);
year_interest = r*12

当然了，我们发现，已知递推关系式，完全可以写出递归函数来计算$x_{360}$ 的数学关系式

求递归函数代码：

function y = loan(x0, a, r, k)
%功能：已知初始欠款数x0，每个月还款数a，月利息r，求k个月以后的欠款y
i=0;
y=x0;
for i=1:k
    y=y*(1+r)-a;
end
end

%function y = loan(x0, a, r, k)
%功能：已知初始欠款数x0，每个月还款数a，月利息r
%求k个月以后的欠款y
x0=25.2;
a=0.1436;
%解方程f(r) = loan(x0, a, r, 360) = 0
f=@(r)loan(x0, a, r, 360);
r0=0.001;
r=fzero(f,r0);
year_interest = r*12

这种方法同样可以求出方程的解，但是要注意的是，当数据量巨大时，递归计算会占用很大的空间与时间，这种情况下，则不如推出的关系式方便。同时线性迭代时要注意迭代关系式是否收敛
本题结果为：

year_interest =
    0.0553

三、斜抛问题

已知一颗炮弹发射做斜抛运动，初速度为200米每秒，击中水平距离360米时，垂直距离160m的目标，问当忽略空气阻力时发射角多大

这是一道简单的高中物理问题，下面我们建立数学模型来解决它

分析
我们对炮弹的轨迹正交分解，设发射角为alpha，用x(t),y(t)表示 t 时刻炮弹的位置，初始速度为v0，由于炮弹发射出去只受到重力，我们可以列出水平与竖直两个方向的方程：
$v_0$*cos(alpha) *t = 360
$v_0$*sin(alpha) *t - $\frac{1}{2}$ g t² = 160
然后，我们就可以利用非线性方程组求解了

代码如下：

%炮弹发射为斜抛运动，已知初速度为200米每秒，
%问击中水平距离360米时，垂直距离160m的目标，问当忽略空气阻力时发射角多大
%设发射角为alpha，用x(t),y(t)表示 t时刻炮弹的位置，初始速度为v0，由于炮弹发射出去只受到重力
%X=[alpha, T]
g=9.8;
fun = @(X) [200*cos(X(1)) * X(2)- 360
    200*sin(X(1)) * X(2) - 1 / 2*g*X(2)^2 - 160];
X0=[pi/6, 5];  %估计初值
X=fsolve(fun,X0)
alpha = X(1) / pi *180 %弧度制转换角度值

当然了，我们也可以动态的观察炮弹的运动轨迹

%炮弹运动轨迹
g=9.8;
f=@(a) 180*sin(2*a)-g/2*1.8^2 - 160*cos(a)^2;
a=fzero(f, pi/6);
figure
axis([0,360,0,160])
hold on
for t=0:0.001:2.012
    x=200*cos(a)*t;
    y=200*sin(a)*t - 1/2*g*t^2;
    plot(x,y,'r.')
    pause(0.005)
end

结果为：

X =
    0.4633    2.0121
alpha =
   26.5444

---

总结

关于数学建模的大致过程基本可归纳为（个人理解），
梳理关系式——写出数学方程或方程组——画出相关函数曲线分析问题——解方程得到问题的解——回带检验误差