Задача о траектории частицы (неправильно работает программа)

@ZergsInPanic · Регистрация: 25.01.2016

Студворк — интернет-сервис помощи студентам

Здравствуйте, речь идет о решении задачи о движении в магнитном поле и в поле силы сопротивления зависящей от скорости, набранный мною код в матлабе считается без ошибок, но считается неверно - выдает не ожидаемый результат) Вот код моих файлов в матлабе:
fun.m:

Matlab M
1
2
3
function underint = fun(t)
underint=exp(-((t).^2));
end

xforce.m:

Matlab M
1
2
3
4
5
6
function eq = xforce(vx,vy)
global alfa beta gama c;
I = quad(@fun,0,beta*sqrt(vx.^2+vy.^2));
eq = (alfa.*vx./(sqrt((vx.^2+vy.^2).^3))).*((2/sqrt(pi)).*I)-c.*sqrt(vx.^2+vy.^2).*exp(-gama.*((vx.^2+vy.^2)));
%написать уравнение задающее силу на каждую компоненту силы для каждого компаненты уравнения.
end

yforce.m:

Matlab M
1
2
3
4
5
6
function eq = yforce(vx,vy)
global alfa beta gama c v0;
I = quad(@fun,0,beta*sqrt(vx.^2+vy.^2));
eq = (alfa.*vy./(sqrt((vx.^2+vy.^2).^3))).*((2/sqrt(pi)).*I)-c.*sqrt(vx.^2+vy.^2).*exp(-gama.*((vx.^2+vy.^2)));
%написать уравнение задающее силу на каждую компоненту силы для каждого компаненты уравнения.
end

diffeq.m:

Matlab M
1
2
3
4
5
6
7
function deq = diffeq(t,r)
deq = zeros(4,1);
deq(1) = r(1);
deq(2) = r(2);
deq(3) = r(2)+xforce(r(1),r(2));
deq(4) = r(1)+yforce(r(1),r(2));
end

result.m:

Matlab M
1
2
3
4
5
6
7
8
9
10
clear all;
global alfa beta gama c v0;
v0 = 0.25;
alfa = 1/38.2;
R0 = 0.1;
gama = 54.4;
beta = sqrt(gama);
c = 2*sqrt(gama/pi);
[T,R] = ode45(@diffeq,[0 50],[0 0 v0 0]);
plot(T,R);

Честно говоря это не мое задание, своя я сделал оно было внятнее, это для одногрупницы так что от нее я думаю вам еще большее спасибо) Если есть у кого какие идеи как получить реальную траекторию частицы движущейся в цилендрическом поле то буду очень благодарен!
P.S.
Только сейчас осознал что сила трения заканчивает свое действие на R0, пойду дорабатывать) но реально результат должен хотя бы начаться красиво, так что помогите мне чем можете) Заранее спасибо!

@ZergsInPanic · 25.01.2016, 12:27 **[ТС]**

Если кому то интересно вот весь текст задания сфотоный с листка:

@ZergsInPanic · 25.01.2016, 12:38 **[ТС]**

Оказывается если посмотреть на вектора T и R выходящие из решения дифф уравнения то скорости не определяются и соответственно траектория не ищется( Что может быть не так?

@SSC · 25.01.2016, 13:06

1. Для deq(3) и deq(4) у Вас странные выражения.
То что записано у Вас примерно так:
Изменение скорости по данному направлению равно скорости по другому направлению плюс некоторая величина, во всяком случае единицы измерения даже не сходятся.
2. В условиях задачи и решении отсутствует масса частицы, а ведь dV/dt=F/m, а в задании dV/dt=...+fт.
3. В исходных данных у Вас есть R0, а в решении диффуравнения эта величина не используется.
4. Как Вы выполняете условие по ограничению величины шага интегрирования dt<=0.2

Не по теме:

В общем подставили Вы одногрупницу ;D

@bobah16 · 25.01.2016, 13:58

ZergsInPanic, интеграл с предынтегральным множителем в силе сопротивления можно заменить на erf(x) (встроенная в матлаб функция).

@ZergsInPanic · 26.01.2016, 09:49 **[ТС]**

Сообщение от SSC

В условиях задачи и решении отсутствует масса частицы, а ведь dV/dt=F/m, а в задании dV/dt=...+fт

Понимаю вас прекрасно, видимо в задании предполагается обезразмеренная сила, да с законной точки зрения действительно странно, но даже так F=dP/dt)

Сообщение от SSC

измерения даже не сходятся.

Тоже понимаю) Но если посмотреть на задание то там даже вектор магнитной напряженности безразмерен и он же в векторном произведении и оставляет компоненту скорости помноженную на единицу)
В любом случае спасибо что наставления и что провели время за распознаванием моего кода)
P.S. А со временем просто, в матлабе точность решения гораздо ниже верхней границы точности решения указанного в задании)

Добавлено через 18 часов 35 минут
Тема все еще актуальна!

Добавлено через 33 секунды
Тема все еще актуальна!

@SSC · 26.01.2016, 11:34

Задание стало не видно, поэтому по памяти:
Модуль diffeq.m:
В него надо добавить - ( а так же в модуль result.m: )

Matlab M
1
global R0

Сообщение от ZergsInPanic

deq(3) = r(2)+xforce(r(1),r(2));
deq(4) = r(1)+yforce(r(1),r(2));

Если считать, что функции xforce, yforce определяют составдяющие тормозящей силы
И если сила действующая от постоянного поля пропорциональна скорости по другой оси (это проверь пишу просто предположение, т.к. у Вас вообще от координаты), тогда вместо этих строк надо поставить

Matlab M
1
2
3
4
5
6
7
if r(1)^2+r(2)^2<=R0^2
  deq(3)=r(4)+xforce(r(1),r(2));
  deq(4)=r(3)+yforce(r(1),r(2));
else
  deq(3)=r(4);
  deq(4)=r(3);
end

@ZergsInPanic · 26.01.2016, 12:04 **[ТС]**

Сообщение от SSC

тогда вместо этих строк надо поставить

Пока еще не пробовал, но можете объяснить каким образом появляются дополнительные компоненты r(3) и r(4) если скорости только vx и vy?

Добавлено через 19 минут

Сообщение от SSC

а так же в модуль result.m:

Что то вновь не физчный результат получился. По логике должна быть лента выглядящая как пружина сжатая , или если математическим языком кривая подобная r(phi)=phi в полярных координатах.

@SSC · 26.01.2016, 14:08

Сообщение от ZergsInPanic

plot(T,R);

Сообщение от ZergsInPanic

По логике должна быть лента выглядящая как пружина сжатая

Какой из параметров R должен быть в зависимости от времени t "как пружина сжатая"?
Скорее всего Вы говорите о траектории движения частицы, то есть зависимости y(x), которая не выводится Вашей командой plot.
Ваше задание я вижу как на картинке внизу, по памяти я проверить ничего не могу.
У Вас перепутаны знаки и номера переменных.
Приведенный ниже код считает при отсутствии силы сопротивления ( R0=0 ). Траектория движения круговая, постепенно, из-за ошибок округления стягивающаяся к центру окружности.
Если R0=0.1, то рсчет идет только до t=0.34, дальше программа виснет из-за ошибок в определении тормозящей силы.
В функциях xforce, уforce I становится маленькой, а величины eq достаточно большими и постоянно меняющими знак.
Поэтому ищи там ошибку.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function reshenie
clear all;clc;close all;
global alfa beta gama c v0;
global R0 
v0 = 0.25;
alfa = 1/38.2;
R0 = 0.1;
R0 = 0.0;
gama = 54.4;
beta = sqrt(gama);
c = 2*sqrt(gama/pi);
%options=odeset('RelTol',1e-5, 'MaxStep', 0.1);
%[T R] = ode45(@diffeq,[0:0.01:0.34],[0 0 v0 0], options);
[T R] = ode45(@diffeq,[0:0.05:36],[0 0 v0 0]);
figure 
plot(T,R(:,1),T,R(:,2),T,R(:,3),T,R(:,4));
figure
plot(R(:,1),R(:,2));
end
 
function deq = diffeq(t,r)
global R0 
% Для r 1-x  2-y  3-Vx  4-Vy
deq = zeros(4,1);
deq(1) = r(3);
deq(2) = r(4);
if r(1)^2+r(2)^2<R0^2
    deq(3) = -r(4)-xforce(r(3),r(4));
    deq(4) = r(3)-yforce(r(3),r(4));
else
    deq(3) = -r(4);
    deq(4) = r(3);
end
end
 
function eq = xforce(vx,vy)
global alfa beta gama c;
I = quad(@fun,0,beta*sqrt(vx.^2+vy.^2));
eq = (alfa.*vx./(sqrt((vx.^2+vy.^2).^3))).*((2/sqrt(pi)).*I)-c.*sqrt(vx.^2+vy.^2).*exp(-gama.*((vx.^2+vy.^2)));
%написать уравнение задающее силу на каждую компоненту силы для каждого компаненты уравнения.
end
 
function eq = yforce(vx,vy)
global alfa beta gama c v0;
I = quad(@fun,0,beta*sqrt(vx.^2+vy.^2));
eq = (alfa.*vy./(sqrt((vx.^2+vy.^2).^3))).*((2/sqrt(pi)).*I)-c.*sqrt(vx.^2+vy.^2).*exp(-gama.*((vx.^2+vy.^2)));
%написать уравнение задающее силу на каждую компоненту силы для каждого компаненты уравнения.
end
 
function underint = fun(t)
underint=exp(-((t).^2));
end

@ZergsInPanic · 28.01.2016, 12:30 **[ТС]**

Сообщение от SSC

У Вас перепутаны знаки и номера переменных.

Спасибо большое очень помогли, надеюсь одногрупнице поставят 5)
Всего вам доброго и удачи)

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:	Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .	Как дизайн сайта влияет на конверсию: 7 решений, которые реально повышают заявки Neotwalker 08.03.2026 Многие до сих пор воспринимают дизайн сайта как “красивую оболочку”. На практике всё иначе: дизайн напрямую влияет на то, оставит человек заявку или уйдёт через несколько секунд. Даже если у вас. . .	Модульная разработка через nuget packages DevAlt 07.03.2026 Сложившийся в . Net-среде способ разработки чаще всего предполагает монорепозиторий в котором находятся все исходники. При создании нового решения, мы просто добавляем нужные проекты и имеем. . .
Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .	SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip Сканируйте QR-код на мобильном и вы увидите, что появится джойстик для управления главным героем. . . .

@ZergsInPanic 0 / 0 / 0 Регистрация: 25.01.2016 Сообщений: 13
	25.01.2016, 12:27 [ТС]
	Если кому то интересно вот весь текст задания сфотоный с листка: 0

@ZergsInPanic 0 / 0 / 0 Регистрация: 25.01.2016 Сообщений: 13
	25.01.2016, 12:38 [ТС]
	Оказывается если посмотреть на вектора T и R выходящие из решения дифф уравнения то скорости не определяются и соответственно траектория не ищется( Что может быть не так? 0

@SSC $Эксперт по математике/физике$ 3390 / 1913 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	25.01.2016, 13:06
	1. Для deq(3) и deq(4) у Вас странные выражения. То что записано у Вас примерно так: Изменение скорости по данному направлению равно скорости по другому направлению плюс некоторая величина, во всяком случае единицы измерения даже не сходятся. 2. В условиях задачи и решении отсутствует масса частицы, а ведь dV/dt=F/m, а в задании dV/dt=...+fт. 3. В исходных данных у Вас есть R0, а в решении диффуравнения эта величина не используется. 4. Как Вы выполняете условие по ограничению величины шага интегрирования dt<=0.2 Не по теме: В общем подставили Вы одногрупницу ;D 1

@bobah16 373 / 343 / 42 Регистрация: 14.07.2015 Сообщений: 2,890
	25.01.2016, 13:58
	ZergsInPanic, интеграл с предынтегральным множителем в силе сопротивления можно заменить на erf(x) (встроенная в матлаб функция). 1