# 拟合常微分方程 (ODE)

### 洛仑兹方程组：定义和数值解

`$\begin{array}{l}\frac{dx}{dt}=\sigma \left(y-x\right)\\ \frac{dy}{dt}=x\left(\rho -z\right)-y\\ \frac{dz}{dt}=xy-\beta z.\end{array}$`

```sigma = 10; beta = 8/3; rho = 28; f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)]; xt0 = [10,20,10]; [tspan,a] = ode45(f,[0 100],xt0); % Runge-Kutta 4th/5th order ODE solver figure plot3(a(:,1),a(:,2),a(:,3)) view([-10.0 -2.0])```

```[tspan,a] = ode45(f,[0 1/10],xt0); % Runge-Kutta 4th/5th order ODE solver figure plot3(a(:,1),a(:,2),a(:,3)) view([-30 -70])```

### 对 ODE 解进行圆周路径拟合

• 路径与 x-y 平面间的夹角 $\theta \left(1\right)$

• 平面与 x 轴向倾角间的夹角 $\theta \left(2\right)$

• 半径 R

• 速度 V

• 相对于时间 0 的时移 t0

• 空间 delta 中的三维位移

`type fitlorenzfn`
```function f = fitlorenzfn(x,xdata) theta = x(1:2); R = x(3); V = x(4); t0 = x(5); delta = x(6:8); f = zeros(length(xdata),3); f(:,3) = R*sin(theta(1))*sin(V*(xdata - t0)) + delta(3); f(:,1) = R*cos(V*(xdata - t0))*cos(theta(2)) ... - R*sin(V*(xdata - t0))*cos(theta(1))*sin(theta(2)) + delta(1); f(:,2) = R*sin(V*(xdata - t0))*cos(theta(1))*cos(theta(2)) ... - R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2); ```

```lb = [-pi/2,-pi,5,-15,-pi,-40,-40,-40]; ub = [pi/2,pi,60,15,pi,40,40,40]; theta0 = [0;0]; R0 = 20; V0 = 1; t0 = 0; delta0 = zeros(3,1); x0 = [theta0;R0;V0;t0;delta0]; [xbest,resnorm,residual] = lsqcurvefit(@fitlorenzfn,x0,tspan,a,lb,ub);```
```Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. ```

```soln = a + residual; hold on plot3(soln(:,1),soln(:,2),soln(:,3),'r') legend('ODE Solution','Circular Arc') hold off```

```figure plot3(a(:,1),a(:,2),a(:,3),'b.','MarkerSize',10) hold on plot3(soln(:,1),soln(:,2),soln(:,3),'rx','MarkerSize',10) legend('ODE Solution','Circular Arc') hold off```

### 对圆弧进行 ODE 拟合

`type paramfun`
```function pos = paramfun(x,tspan) sigma = x(1); beta = x(2); rho = x(3); xt0 = x(4:6); f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)]; [~,pos] = ode45(f,tspan,xt0); ```

```xt0 = zeros(1,6); xt0(1) = sigma; xt0(2) = beta; xt0(3) = rho; xt0(4:6) = soln(1,:); [pbest,presnorm,presidual,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,soln);```
```Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. ```

`fprintf('New parameters: %f, %f, %f',pbest(1:3))`
```New parameters: 9.132446, 2.854998, 27.937986 ```
`fprintf('Original parameters: %f, %f, %f',[sigma,beta,rho])`
```Original parameters: 10.000000, 2.666667, 28.000000 ```

```figure hold on odesl = presidual + soln; plot3(odesl(:,1),odesl(:,2),odesl(:,3),'b') plot3(soln(:,1),soln(:,2),soln(:,3),'r') legend('ODE Solution','Circular Arc') view([-30 -70]) hold off```

### ODE 拟合中的问题

```options = optimoptions('lsqcurvefit','FiniteDifferenceStepSize',1e-4,... 'FiniteDifferenceType','central'); [pbest2,presnorm2,presidual2,exitflag2,output2] = ... lsqcurvefit(@paramfun,xt0,tspan,soln,[],[],options);```
```Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. ```

`disp([presnorm,presnorm2])`
``` 20.0637 20.0637 ```