问题描述
我正在尝试使用Runge-Kutta四阶方法对摆进行简单的模拟。我正在使用p5.js。通常，它可以正确地计算角度，但有时它只是开始随机旋转等。我不知道如何解决这个问题，我认为我说的算法的实现似乎有问题吗？
我以陶韬的《计算物理入门》中的代码为灵感，看起来很相似
上述代码在Fortran中的第一部分
第二部分
let screenWidth = 1300;
let screenHight = 1970;
let angleChangeDifference;
let gSlider;
let lSlider;
let aFrequencySlider;
let dumpingSlider;
let startAngleSlider;
let timestepSlider;
let timeMaxSlider;
let dForceSlider;
let initForceSlider;
let startAngle = 0;
let currentAngle = 0;
let circleX = 0;
let circleY = 200;
let circleRWidth = 100;
let circleRHeight = 100;
let lineXStart = 0;
let lineYStart = 0;
let lineXEnd = circleX;
let lineYEnd = (circleY - circleRHeight / 2);
function setup() {
createCanvas(screenWidth,screenHight);
changeScreenDeafultStartingPoint(screenWidth / 2,100);
frameRate(60)
createSliders();
setInterval(showSliderValue,100);
startButton.mouseClicked(start);
restartButton.mouseClicked(restart);
chartButton.mouseClicked(enableChart);
background(200);
}
let angleSign = '\u00B0';
let omegaSign = '\u03C9';
let chartOn = false;
let step = 0;
function draw() {
startupConfiguration()
showSliderValue()
if (step >= 1) {
rotatePendulum();
if (step == 1) {}
if (step == 2)
startButton.remove();
chartButton.position(20,340);
if (chartOn == true) {
createChart(0,1,'czas (s)','kat (' + angleSign + ')',degreesArr);
createChart(0,320,omegaSign + ' (' + angleSign + '/s)',omegaArr);
scale(2);
createPhaseChart(290,80,degreesArr);
scale(0.5);
}
}
line(lineXStart,lineYStart,lineXEnd,lineYEnd + (20 * lSlider.value()));
fill(200,76,43)
ellipse(circleX,circleY + (20 * lSlider.value()),circleRWidth,circleRHeight);
}
function createSliders() {
gSlider = createSlider(0.05,20,9.81,0.01);
gSlider.position(1100,-90);
lSlider = createSlider(0.5,10,5,0.5);
lSlider.position(1100,-50);
initForceSlider = createSlider(-5,0.5);
initForceSlider.position(1100,50);
dForceSlider = createSlider(-2,2,0.9,0.05)
dForceSlider.position(1100,90);
aFrequencySlider = createSlider(-2,2 / 3,1 / 3);
aFrequencySlider.position(1100,130);
dumpingSlider = createSlider(0.00,1.5,0.5,0.05);
dumpingSlider.position(1100,170);
startAngleSlider = createSlider(-Math.PI / 2,Math.PI / 2,Math.PI / 32);
startAngleSlider.position(1100,210);
timestepSlider = createSlider(0,1000,100,10);
timestepSlider.position(1100,250);
timeMaxSlider = createSlider(10,10000,10);
timeMaxSlider.position(1100,290);
startButton = createButton('ZATWIERDZ',false);
startButton.position(100,310)
restartButton = createButton('RESTART',false);
restartButton.position(20,310)
chartButton = createButton('WYKRES',false);
chartButton.position(-200,340);
}
function showSliderValue() {
background(200);
fill(0,0)
text('sila poczatkowa',440,-60)
text(initForceSlider.value(),400,-42)
text('sila sprawcza',-20)
text(dForceSlider.value(),-2)
text('czestosc katowa',20)
text(aFrequencySlider.value(),42)
text('tlumienie',60)
text(dumpingSlider.value(),82)
text('kat poczatkowy',100)
text(int(degrees(startAngleSlider.value())),122)
text('krok czasowy (N1)',140)
text(timestepSlider.value(),162)
text('dlugosc symulacji (N2)',180)
text(timeMaxSlider.value(),202)
}
function start() {
angleIndex = 0;
step++;
startAngle = startAngleSlider.value();
currentAngle = startAngle;
angleChangeDifference = simulate();
rotatePendulum()
startButton.html("START")
}
function restart() {
window.location.reload();
}
function enableChart() {
chartOn = true;
}
function createChart(moveByX,moveByY,xName,yName,table) {
rotate(-currentAngle);
scale(1.1);
translate(moveByX,moveByY);
strokeWeight(1);
line(-500,500,530,500);
line(-500,700,700);
line(-500,-500,700);
line(530,700);
strokeWeight(1);
let counter = 0;
for (i = 510; i < 700; i += 10) {
if (counter < 9 && counter % 2 == 0)
text(90 - 10 * counter,-520,i + 5)
else if (counter == 9 && counter % 2 == 0)
text(90 - 10 * counter,-515,i + 5)
else if (counter > 9 && counter % 2 == 0)
text(90 - 10 * counter,-525,i + 5)
line(-505,i,i);
counter++;
}
textSize(25);
text(xName,-20,750)
textSize(12);
counter = 0;
for (i = -490; i < 535; i += 25) {
line(i,705);
if (counter % 4 == 0) {
line(i,705);
text(counter * 2.5,i - 5,715);
}
counter++;
}
rotate(-90);
textSize(25);
text(yName,-670,-550)
textSize(12);
rotate(90);
fillChartByTableValues(table);
translate(-moveByX,-moveByY);
scale(0.91);
rotate(currentAngle);
}
function fillChartByTableValues(table) {
strokeWeight(2);
stroke(0,255);
for (i = 0; i < timeArr.length - 1; i++) {
FirstPointX = -490 + timeArr[i] * 10;
FirstPointY = 600 + table[i] * (-1);
SecondPointX = -490 + timeArr[i + 1] * 10;
SecondPointY = 600 + table[i + 1] * (-1);
line(FirstPointX,FirstPointY,SecondPointX,SecondPointY);
}
stroke(0,0);
strokeWeight(0.1);
}
function createPhaseChart(moveByX,-300,700);
line(-300,700);
strokeWeight(1);
let counter = 0;
textSize(8);
for (i = 510; i < 700; i += 10) {
if (counter < 9 && counter % 2 == 0)
text(90 - 10 * counter,-517,i + 3)
else if (counter == 9 && counter % 2 == 0)
text(90 - 10 * counter,-512,i + 3)
else if (counter > 9 && counter % 2 == 0)
text(90 - 10 * counter,i + 3)
line(-505,i);
counter++;
}
textSize(12);
textSize(15);
text(xName,-430,735)
textSize(12);
counter = 0;
textSize(8);
for (i = -490; i < -300; i += 10) {
line(i,705);
if (counter < 9 && counter % 2 == 0)
text(-90 + 10 * counter,i - 7,715);
else if (counter == 9 && counter % 2 == 0)
text(-90 + 10 * counter,i - 2,715);
else if (counter > 9 && counter % 2 == 0)
text(-90 + 10 * counter,i - 4,715);
counter++;
}
textSize(12);
rotate(-90);
textSize(15);
text(yName,-620,-528)
textSize(12);
rotate(90);
fillPhaseChartByTableValues(degreesArr,omegaArr)
translate(-moveByX,-moveByY);
scale(0.91);
rotate(currentAngle);
}
function fillPhaseChartByTableValues(tableX,tableY) {
translate(-400,600);
strokeWeight(1);
stroke(0,255);
for (i = 0; i < tableX.length; i++) {
ellipse(tableX[i],tableY[i],0.5);
}
translate(400,-600);
stroke(0,0);
strokeWeight(0.1);
}
function startupConfiguration() {
background(200);
angleMode(degrees);
changeScreenDeafultStartingPoint(screenWidth / 2,100);
}
function changeScreenDeafultStartingPoint(x,y) {
translate(x,y);
}
let angleIndex = 0;
function rotatePendulum() {
currentAngle = angleChangeDifference[angleIndex] * (180 / PI);
rotate(currentAngle);
if (step > 1) {
angleIndex++
}
}
function calculateIntegral(t,q,dt,f) {
let k1 = f(t,q).map(val => val * dt);
let temp = k1.map(val => val * 0.5);
temp = temp.map((val,index) => val + q[index])
let k2 = f(t + 0.5 * dt,temp).map(val => val * dt);
temp = k2.map(val => val * 0.5);
temp = temp.map((val,index) => val + q[index])
let k3 = f(t + 0.5 * dt,temp).map(val => val * dt);
temp = q.map((val,index) => val + k3[index]);
let k4 = f(t + dt,temp).map(val => val * dt);
temp = k2.map((val,index) => val + k3[index])
temp = temp.map(val => val * 2)
temp = temp.map((val,index) => (val + k1[index] + k4[index]) / 6)
temp = temp.map((val,index) => val + q[index])
return [t + dt,temp];
}
function modelPendulum(t,q) {
let c = dumpingSlider.value();
let fw = dForceSlider.value();
let w = aFrequencySlider.value();
let x1 = q[0];
let x2 = q[1];
return [x2,-(Math.sin(x1)) - (c * x2) + (fw * Math.cos(w * t))];
}
let degreesArr,timeArr;
function simulate() {
let t = 0.0;
let dt = (3 * Math.PI) / timestepSlider.value()
let tf = timeMaxSlider.value() * dt
let q = [startAngle,initForceSlider.value()];
let Nt = int(Math.round((tf - t) / dt)) + 1;
let solution = new Array(q.length + 1);
for (i = 0; i < q.length + 1; i++) {
solution[i] = new Array(Nt).fill(0);
}
solution[0][0] = t;
solution[1][0] = q[0];
solution[2][0] = q[1];
k = 1;
while (t <= tf) {
let temporaryResult = [];
temporaryResult = calculateIntegral(t,modelPendulum);
t = temporaryResult[0];
q = temporaryResult[1];
solution[0][k] = t;
solution[1][k] = q[0];
solution[2][k] = q[1];
k = k + 1
}
timeArr = solution[0];
degreesArr = solution[1];
omegaArr = solution[2];
let counter = 0;
let ifChaos = false;
while (counter != degreesArr.length - 1 && ifChaos != true) {
if (degreesArr[counter] > 13.5 || degreesArr[counter] < -13.5) {
ifChaos = true;
}
counter++;
}
if (ifChaos == true) {
degreesArr = degreesArr.map(val => val * 5.32);
omegaArr = omegaArr.map(val => val * 23.32);
} else {
degreesArr = degreesArr.map(val => val * 35.32);
omegaArr = omegaArr.map(val => val * 35.32);
}
return solution[1]
}<script src="https://cdn.jsdelivr.net/npm/p5@1.1.9/lib/p5.js"></script>解决方法
我不熟悉P5，因此这是一个利用ThreeJS利用Mathematics Stack Exchange中的Runge-Kutta算法的解决方案。
为方便起见，我将Runge-Kutta算法包装在一个类中，构造函数采用了以下初始参数：
- 重力加速度常数
g（对于地球，为9.81米/秒/秒）， - 钟摆长度（以米为单位）
- 钟摆的初始角度（0垂直向下）
- 初始角速度（以米/秒为单位），并且
- 最大时间增量。 （由于使用了Runge-Kutta来求解以时间为变量的二阶微分方程，因此根据作者的实验，人们似乎无法过度延长增量时间增量，而仍然保持所产生的摆位置和速度的精度。此参数只需限制在
t方法中传递的最大updatePosition值，默认值为0.1s。）
为帮助使用Runta-Kutta算法，下面的代码模拟了两个1米的摆锤：
- 第一个具有-PI / 2（-90度）初始位置且没有起始角速度的机器人。
- 第二个具有初始位置PI（180度），起始角速度很小。
<script type="module">
import * as THREE from 'https://cdn.jsdelivr.net/npm/three@0.115.0/build/three.module.js';
class RungeKutta {
constructor( g,pendulumLength,initialAngle,angularVelocity,maxTimeDelta ) {
this.g = g;
this.pendulumLength = pendulumLength;
this.theta = initialAngle;
this.omega = angularVelocity;
this.maxTimeDelta = maxTimeDelta || 0.1;
}
updatePosition( t ) {
let self = this;
function omegaDot( theta ){
return -( self.g / self.pendulumLength ) * Math.sin( theta );
}
function thetaDot( omega ){
return omega;
}
// If the browser tab becomes inactive,then there will be a large
// time delta,which will disrupt the RungeKutta algorithm. If more
// than max allowed seconds has lapsed,then reset the timer.
if ( self.maxTimeDelta < t ) {
t = self.maxTimeDelta;
}
let aomega = omegaDot( self.theta );
let atheta = thetaDot( self.omega );
let bomega = omegaDot( self.theta + 0.5 * t * atheta );
let btheta = thetaDot( self.omega + 0.5 * t * aomega );
let comega = omegaDot( self.theta + 0.5 * t * btheta );
let ctheta = thetaDot( self.omega + 0.5 * t * bomega );
let domega = omegaDot( self.theta + t * ctheta );
let dtheta = thetaDot( self.omega + t * comega );
self.omega = self.omega + ( t / 6 ) * ( aomega + 2 * bomega + 2 * comega + domega );
self.theta = self.theta + ( t / 6 ) * ( atheta + 2 * btheta + 2 * ctheta + dtheta );
return self;
}
}
//
// Set up the ThreeJS environment.
//
var renderer = new THREE.WebGLRenderer();
renderer.setSize( window.innerWidth,window.innerHeight );
document.body.appendChild( renderer.domElement );
var camera = new THREE.PerspectiveCamera( 45,window.innerWidth / window.innerHeight,1,500 );
camera.position.set( 0,100 );
camera.lookAt( 0,0 );
var scene = new THREE.Scene();
//
// Create the pendulum mesh.
//
var length = 30,width = 1;
var shape = new THREE.Shape();
shape.moveTo( -width / 2,0 );
shape.lineTo( +width / 2,-length );
shape.lineTo( -width / 2,0 );
var extrudeSettings = {
steps: 2,depth: 2,bevelEnabled: true,bevelThickness: .25,bevelSize: .25,bevelOffset: 0,bevelSegments: 1
};
var geometry = new THREE.ExtrudeBufferGeometry( shape,extrudeSettings );
var material = new THREE.MeshBasicMaterial( { color: 0x00ff00 } );
var mesh0 = new THREE.Mesh( geometry,material );
mesh0.position.x = -17;
scene.add( mesh0 );
var mesh1 = mesh0.clone();
mesh1.position.x = +17;
scene.add( mesh1 );
//
// And now animate the pendulum using RungeKutta.
//
let pendulumState0 = new RungeKutta( 9.81,-Math.PI / 2,0.1 );
let pendulumState1 = new RungeKutta( 9.81,Math.PI,0.01,0.1 );
let now = performance.now();
let lastTimer = now;
var animate = function () {
requestAnimationFrame( animate );
now = performance.now();
pendulumState0.updatePosition( ( now - lastTimer ) / 1000 );
pendulumState1.updatePosition( ( now - lastTimer ) / 1000 );
lastTimer = now;
mesh0.rotation.z = pendulumState0.theta;
mesh1.rotation.z = pendulumState1.theta;
renderer.render( scene,camera );
};
animate();
</script>
希望这将有助于您实施P5。
,这是pendulum与{5}一起使用Leapfrog method的示例，但是可以很容易地适应使用Runge-Kutta 4（Leapfrog具有辛的优势，与RK有所不同）。
对于一个简单的摆，您必须求解二阶微分方程（牛顿或欧拉-拉格朗日方程）
d ^ 2 O / dt ^ 2 =-（g / l）sin（O）。
其中O是与垂直线的夹角。
首先应用数值方法很方便，可以将其转换为两个一阶微分方程的等效系统。
dw / dt =-（g / l）sin（O）
dO / dt = w
例如，您可以使用Runge-Kutta4或其他方法来解决此系统。一个非常简单的方法就是简单的二阶Leapfrog方法。您将时间离散化为dt这样
w_n == w（（n-1 / 2）dt）
O_n == O（n dt）
，然后循环循环执行以下操作
w_ {n + 1} = w_n-（g / l）sin（O_n）dt
O_ {n + 1} = O_n + w_ {n + 1} dt
此方法在此code的钟摆类中实现，及其参数，p5中用于可视化的非常简单的render（）函数以及例如用于计算能量的方法（可用于检查节约能源（如果不强制）。这非常简单，您可以对其进行测试。当然，应使用较小的时间步长dt以便获得较高的精度。
如果您想使用Runge-Kutta4或其他高阶方法，可以方便地编写函数Fw（t，w，O）和FO（t，w，O）以便方程式读取
dw / dt = Fw（t，w，O）（==-（g / l）sin（O）+ F0 cos（W t））
dO / dt = FO（t，w，O）（== w）
现在我在其中添加了一个强制F0 cos（W t），以示以防万一，如何添加任何与时间有关的力。
然后您可以及时离散化（没有越级跨步）
t_n == n dt
w_n == w（n dt）
O_n == O（n dt）
并在每个步骤n中计算RK4的数量k1，k2，k3，k4。 请注意，在这种情况下，ki是二维向量{kiw，kiO}，函数f（t，w，O）== {Fw（t，w，O），FO（t，w，O）}为也是二维向量。
考虑到这一点，并分别定义函数Fw和FO，您可以用RK4方法轻松替换p5示例1的Leapfrog方法，并根据需要添加强制。
如果您看到非常奇怪的行为，则应检查标志和RK4方法的正确设置。
我希望这会有所帮助。
let P1;
function setup() {
createCanvas(720,400);
P1=new PenduloSimple(3.14,0.0);
}
function draw() {
background(220);
P1.render();
P1.leapFrog();
}
class PenduloSimple {
constructor(ang,velang) {
this.ang = ang;
this.velang = velang;
this.g = 9.8;
this.dt = 0.01;
this.l = 1.0;
this.m = 0.1;
this.E0 = this.Energy();
}
Energy() {
let E=
this.m * this.l * this.l * this.velang * this.velang * 0.5 - this.m * this.g * this.l * cos(this.ang);
return E;
}
leapFrog() {
// Método de Leapfrog
this.velang =
this.velang + this.dt*(-this.g/this.l)*sin(this.ang);
this.ang = this.ang + this.dt*this.velang;
}
render(){
var x0 = width / 2;
var y0 = height / 2;
var mult0 = 100;
var x = this.l * sin(this.ang);
var y = this.l * cos(this.ang);
var xplot = x0 + x * mult0;
var yplot = y0 + y * mult0;
// Draw a circle
stroke(50);
fill(100);
ellipse(xplot,yplot,24,24);
// La cuerda
stroke(50);
line(x0,y0,xplot,yplot);
let E = this.Energy();
stroke(0);
fill(255,0);
rect(10,height * 0.5,20,-E * height * 0.45 / this.E0);
fill(0,255,0);
rect(30,-this.ang * height * 0.4 / 6.28);
fill(0,255);
rect(50,-this.velang * mult0 * 0.1);
fill(1);
stroke(255);
fill(255,0);
text("E: " + E,width * 0.1,20);
fill(50,200,50);
text("ang: " + this.ang,30);
fill(0,255);
text("velang: " + this.velang,40);
}
}