常微分方程
1、四阶经典步长Runge-Kutta法求解高阶常微分方程
clear;
clc;
% 四阶经典定步长Runge-Kutta法求解高阶常微分方程组
% function
f = @(x,y,z)z;
g = @(x,y,z)2*(y^3);
% 条件
h = 0.01;
len = 0.5/h;
x(1) = 1;
y(1) = -1;
z(1) = -1;
i = 1;
k = 0;
while k < len
K1 = f(x(i),y(i),z(i));
L1 = g(x(i),y(i),z(i));
K2 = f(x(i)+h/2,y(i)+h*K1/2,z(i)+h*L1/2);
L2 = g(x(i)+h/2,y(i)+h*K1/2,z(i)+h*L1/2);
K3 = f(x(i)+h/2,y(i)+h*K2/2, z(i)+h*L2/2);
L3 = g(x(i)+h/2,y(i)+h*K2/2, z(i)+h*L2/2);
K4 = f(x(i)+h,y(i)+h*K3,z(i)+h*L3);
L4 = g(x(i)+h,y(i)+h*K3,z(i)+h*L3);
y(i+1) = y(i)+(K1+2*K2+2*K3+K4)*h/6;
z(i+1) = z(i)+(L1+2*L2+2*L3+L4)*h/6;
x(i+1) = x(i) + h;
i = i + 1;
k = k + 1;
end
plot(x, y, '-o');
title('第一题图');
y = y';
z = z';
2、用ode15s函数、ode23s函数和ode23tb函数求刚性方程并画图比较：
https://ww2.mathworks.cn/help/matlab/ref/ode15s.html
function dydt = vdp1000(t,y)
%VDP1000 Evaluate the van der Pol ODEs for mu = 1000.
%
% See also ODE15S, ODE23S, ODE23T, ODE23TB.
% Jacek Kierzenka and Lawrence F. Shampine
% Copyright 1984-2014 The MathWorks, Inc.
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];
使用 ode15s 求解器对刚性方程组求解，然后绘制解 y 的第一列对时间点 t 的图。ode15s 求解器通过刚性区域所需的步数远远少于 ode45。
%% ode15s
[t1,y1] = ode15s(@vdp1000,[0 3000],[2 0]);
figure(1)
plot(t1,y1(:,1),'-o')
title('ode15s作图')
%% ode23s
[t2,y2] = ode23s(@vdp1000,[0 3000],[2 0]);
figure(2)
plot(t2,y2(:,1),'-*')
title('ode23s作图')
%% ode23tb
[t3,y3] = ode23tb(@vdp1000,[0 3000],[2 0]);
figure(3)
plot(t3,y3(:,1),'-.')
title('ode23tb作图')
ode15s.jpg
ode23s.jpg
ode23tb.jpg
3、求解微分方程组
第三题图
y0 = [0, 0, 0.68, 1, -0.5];
[t, y] = ode45(@fun8_3, [-10, 10], y0);
% plot(y(:,4),y(:,1));
plot (t,y);
title('第三题');
% 微分方程组
function f = fun8_3(t,y)
f = [y(2);
y(3);
-3*y(1)*y(3)+2*(y(2))^2 - y(4);
y(5);
-2.1*y(1)*y(5);];
end
-
已知Appolo卫星的运动轨迹（x，y）满足下面的方程：
方程组.jpg
其中条件1.jpg
试在初值条件为
初值条件.jpg的情况下求解，并绘制卫星轨迹图。
解：
函数appollo，初始化方程组，设定参数
function dx = appollo(t,x)
% mui = μ
% mustar = λ
mui = 1/82;
mustar = 1 - mui;
% 条件r1,r2
% x(1) = x;
% x(3) = y;
r1 = sqrt((x(1) + mui)^2 + x(3)^2);
r2 = sqrt((x(1) - mustar)^2 + x(3)^2);
% 方程组
% dy/dt = x(4)
% d^2y/dt^2 = x(2)
dx = [x(2)
2*x(4) + x(1) - mustar*(x(1)+mui)/r1^3 - mui*(x(1) - mustar)/r2^3
x(4)
-2*x(2) + x(3) - mustar*x(3)/r1^3 - mui*x(3)/r2^3];
主函数
clear;
clc;
% test 4
% x0(i)对应与xi的初值
% x(0) = 1.2;
% x'(0) = 0;
% y(0) = 0;
% y'(0) = -1.04935371
x0=[1.2;0;0;-1.04935371];
options=odeset('reltol',1e-8);
% 该命令的另一种写法是options=odeset;options.reltol=1e-8;
% 是解微分方程时的选项设置，'RelTol',1e-4,是相对误差设置
% 'AbsTol',[1e-4 1e-4 1e-5]是绝对误差设置
tic
[t,y]=ode45(@appollo,[0,20],x0,options);
%t是时间点，y的第i列对应xi的值，t和y的行数相同
toc
plot(y(:,1),y(:,3))
%绘制x1和x3，也就是x和y的图形
title('Appollo卫星运动轨迹')
xlabel('X')
ylabel('Y')
运行结果：
>> test4
时间已过 0.039310 秒。
appollo.jpg
偏微分方程求解工具的学习
图片1.jpg
图片2.jpg
图片3.jpg
图片4.jpg
图片5.jpg
图片6.jpg
图片7.png
图片8.png
