问题描述
我正在尝试使用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/[email protected]/lib/p5.js"></script>解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）