Методика решения физических задач в среде МатЛаб

Введение
1. Среда MATLAB
2. Основные этапы для решения физических задач в MATLAB
3. Физические задачи в MATLAB
Применение численных методов. Примеры готовых задач
4. Решение физических задач на уроках
5. Задачи для самостоятельной работы
Заключение
Помощь в работе с MATLAB
Литература

Методика решения физических задач в среде МатЛаб

n=3;q=0,1
График №6
n=3;q=0,4
График №7
n=3;q=0,5
График №8
Методы Рунге-Кутты для решения физических задач в MATLAB
Методы Рунге-Кутты в MATLAB реализуют функции ode23 (точность ) и ode45 (точность ). В этих функциях шаг интегрирования h выбирается автоматически так, чтобы гарантировать заданную относительную точность решения (как она задается – см. HELP для ode23). По умолчанию она равна 1е-3.
Стандартная процедура решения ode23. Обращение имеет вид:
[T у1]=ode23(‘fun’, [t0 tmax], y0). Здесь t0 – начальное t, y0 – начальное условие, fun – имя функции, вычисляющей правую часть ДУ, у1 – полученное решение, T – значения аргумента, для которых оно получено. Вместо [t0 tmax] можно задать те значения t, для которых мы хотим получить решения. Например, [t0:h:tmax.]
Зададим функцию правой части нашего ДУ. Это делается для всех задач, которые будут рассмотрены ниже.
Function z=mytest1(T, y)
Z=sqrt(T.^2+y.^2)
Решение физических задач
1.Падение тела в одномерном гравитационном поле с высоты Н, без начальной скорости:
.
Эквивалентная система:
.
Где y(1) – координата, y(2) – скорость. Это массивы типа y(:,1), y(:,2).
Текст программы:
function z=mytgrav(T, y)
z=[ y(2) -1./ y(1)^2]’; – столбец правых частей с k=m
tic; [T y1]=ode23(‘mytgrav’, [0 1.1], [1 0]’); toc
figure(1); subplot(2,1,1); plot(T,y1(:,1),’.-b’); title(‘Высота’)
subplot(2,1,2); plot(T,y1(:,2),’r’); title(‘Скорость’); xlabel(‘t’)
2. Задача о собачьей кривой.
Собака со скоростью w бежит по плоскости (х, у) к хозяину, идущему со скоростью v вдоль оси у. Найти траекторию движения собаки y(x). В начальный момент хозяин находится в точке (0, 0), собака – в точке (1, 0).
ДУ для траектории:
.
Переходим к системе:
.
Решим для v/w=0.5
Текст программы:
function z=mydog(x, y)
z=[y(2); 0.5*sqrt(1+y(2)^2)/x]’; – столбец правых частей с k=m
tic; [x y1]=ode23(‘mydog’, [1 0.001], [0 0]’); toc
figure(3); plot(x,y1(:,1),’.-b’); title(‘Собачья кривая’)
xlabel(‘x’); xlabel(‘y(x)’).
4. Решение физических задач на уроках
1. Биения. Функция биения:
.
Несмотря на простоту формулировки, эта задача очень содержательна.
В процессе исполнения программы на экран выводится участок кривой для интервалов времени .
Текст программы:
a1=1.0; % Амплитуды гармонических
a2=1.0; % колебаний
w1=1.0; % Частоты гармонических
w2=1.2; % колебаний
t0=0; % Начальный момент времени
tm=20; % Конечный момент времени
N=600; % Число точек вывода/расчета
T=tm-t0; % Время вывода биений
dt=T/N; % Шаг по времени
t=t0:dt:tm; % Вектор времени
y=a1*cos(w1*t)+a2*cos(w2*t); % Функция биений
plot(t,y); % Вывод графика
После вывода результата расчета на экран можно изменить масштаб осей с повторным выводом соответствующего участка кривой. Так, в приведенном выше примере, переменная t на графике будет изменяться от 0 до 20. Если мы хотим рассмотреть подробности графика в другом диапазоне (например, по t от 1 до 2), то необходимо ввести в командном окне команду axis ([xmin xmax ymin ymax]), где xmin xmax -диапазон вывода по оси x, а ymin ymax -диапазон вывода по оси y. В результате выполнения этого оператора график будет перестроен в указанном масштабе.
Задание. Другого вида биения можно наблюдать при . Для того чтобы увидеть, как действительно выглядят колебания в этом случае, желательно увеличить время вывода биений (т.е. диапазон изменения переменной t) примерно в 10-20 раз.
2. Рассмотрим колебания маятника:
.
Исследуем точки покоя маятника.
Текст программы:
Правые части: для А=1
function z=mypend(t, y)
z=[ y(2); -sin(y(1)];
Нарисуем фазовый портрет системы (для различных начальных условий):
и график потенциальной энергии:
tmax=18;
figure(1); clf; subplot(2,1,1); hold on;
for k=1:12; y0=0.4*(k-6); % разные начальные условия
[t y]=ode23(‘mypend’, [0 tmax], [2*pi y0]’);
plot(y(:,1), y(:,2)); end % Кривые около точки покоя (2*pi, 0)
xlabel(‘x(t)’); ylabel(‘dy/dt’); title(‘Фазовый портрет маятника’);
hold off
x=[-5:0.1:35]; v=1-cos(x);
subplot(2,1,2); plot(x,v,’r’);
title(‘Потенциальная энергия маятника’)
xlabel(‘x’); ylabel(‘U(x)’);
Вопросы учащимся:
- что можно сказать об устойчивости точек покоя ;
- объясните характер движения маятника для каждой фазовой кривой.
3. Задача об изгибании консольной балки.
Из теории упругости известно ДУ для прогиба y(x) и начальные условия:
Зададим в безразмерных единицах:
.
Пусть момент инерции сечения нарисованной балки от x:
.
Текст программы.
Правые части системы ДУ
function z=mybar(x, y)
l=100; J=5*(1+4*exp(-6*x/1));
J1=-120*exp(-6*x/1)/1;J2=720*exp(-6*x/1)/1/1;
Z=[y(2);y(3);y(4);-2*J1/J*y(4)-J2/J*y(3)|];
Зададим начальные условия
y0=[0 0 0.667e-4 2.5335e-8]’;
Решим
[x y]=ode23(‘mybar’, [0 100], y0);
figure(1); subplot(2,1,1); plot(x,y(:,1)); title(‘Прогиб балки’)