Sympy second order ode

2024/9/8 10:48:51

I have a homogeneous solution to a simple second-order ODE, which when I try to solve for initial values using Sympy, returns the same solution. It should substitute for y(0) and y'(0) and yield a solution without constants, but does not. Here is the code to set up the equation (it is a spring balance equation, k = spring constant and m = mass). Some redundant symbols which I use elsewhere, sorry.

%matplotlib inline
from sympy import *
m,k, x,t,s,T, omega,A,B = symbols('m k x t s T omega A B',float=True)
a = symbols('a', positive=True)
f,y,g,H,delta,S=symbols('f y g H delta S',cls=Function)
Eq1 = Eq(m*diff(y(t),t,2)+k*y(t))
Eq1

The result is (correctly): $y{\left (t \right )} = C_{1} e^{- t \sqrt{- \frac{k}{m}}} + C_{2} e^{t \sqrt{- \frac{k}{m}}}$

y(t)=C1e^(−t√(−k/m))+C2e^(t√(−km)), which also has y_n = c1.cos(√(−k/m)t)+c2.sin(√(−k/m)t).

When this equation is solved analytically, and converting to a solution using sines and cosines with omega = sqrt(-k/m), then c1 = y(0) and c2 = y'(0)/omega. So, while the solution is partially involving complex numbers, of course, dsolve simply returns the original homogeneous equation as above. My code to evaluate the ODE at y(0) and y'(0) is:

Eq1_soln_IVP =dsolve(Eq1,y(t),x0=0, ics={y(0): a, y(t).diff(t).subs(t, 0): a})

I appreciate that dsolve simply may not be able to handle this IVP analytically, but I would be surprised if this is so based on its other capacity. Any help as to how this problem and therefore other analytic second order problems can be solved will be much appreciated. The nub of the question is around the :

ics={y(0): a, y(t).diff(t).subs(t, 0): a}

So the solution I tried, which Dietrich confirms, was :

 #Create IVP for y(0)expr = Eq(Eq1_soln_IVP.rhs.subs(sqrt(-k/m),I*omega),y(0))#Create IVP for y'(0)expr2 = Eq(diff(y(t),t).subs(t,0),expr.lhs.diff(t))#Maps all free variables and solves for each where t = 0.solve([expr.subs(t,0),expr2.subs(t,0)])

While it is "a" solution, this seems a very convoluted way of finding y(t) = y(0)cos(omega*t - phi)...which answers the implicit question about some limitations of this solver and the direct question about how the ics arg is resolved. Thanks.

Answer

The Parameter ics in dsolve() does not really work (Issue 4720), so you have to do the substitutions manually. You could try:

from IPython.display import display
import sympy as sysy.init_printing()  # LaTeX-like pretty printing for IPythont = sy.Symbol("t", real=True)
m, k = sy.symbols('m k', real=True)  # gives C_1 Exp() + C_2 Exp() solution
# m, k = sy.symbols('m k', positive=True)  # gives C_1 sin() + C_2 cos() sol.
a0, b0 = sy.symbols('a0, b0', real=True)
y = sy.Function('y')Eq1 = sy.Eq(m*sy.diff(y(t), t, 2) + k*y(t))
print("ODE:")
display(Eq1)print("Generic solution:")
y_sl0 = sy.dsolve(Eq1, y(t)).rhs  # take only right hand side
display(sy.Eq(y(t), y_sl0))# Initial conditions:
cnd0 = sy.Eq(y_sl0.subs(t, 0), a0)  # y(0) = a0
cnd1 = sy.Eq(y_sl0.diff(t).subs(t, 0), b0)  # y'(0) = b0#  Solve for C1, C2:
C1, C2 = sy.symbols("C1, C2")  # generic constants
C1C2_sl = sy.solve([cnd0, cnd1], (C1, C2))# Substitute back into solution:
y_sl1 = sy.simplify(y_sl0.subs(C1C2_sl))
print("Solution with initial conditions:")
display(sy.Eq(y(t), y_sl1))
https://en.xdnf.cn/q/72691.html

Related Q&A

Bulk update using Peewee library

Im trying to update many records inside a table using Peewee library. Inside a for loop, i fetch a single record and then I update it but this sounds awful in terms of performance so I need to do the u…

Can you specify variance in a Python type annotation?

Can you spot the error in the code below? Mypy cant.from typing import Dict, Anydef add_items(d: Dict[str, Any]) -> None:d[foo] = 5d: Dict[str, str] = {} add_items(d)for key, value in d.items():pr…

Django loaddata error

I created a "fixtures" folder in the app directory and put data1.json in there.This is what is in the file:[{"firm_url": "http://www.graychase.com/kadam", "firm_name&…

Parsing JSON string/object in Python

Ive recently started working with JSON in python. Now Im passing a JSON string to Python(Django) through a post request. Now I want to parse/iterate of that data. But I cant find a elegant way to parse…

Removing NaNs in numpy arrays

I have two numpy arrays that contains NaNs:A = np.array([np.nan, 2, np.nan, 3, 4]) B = np.array([ 1 , 2, 3 , 4, np.nan])are there any smart way using numpy to remove the NaNs in b…

Run Python + OpenCV + dlib in Azure Functions

I have created an image processing script in Python (with dlib and OpenCV) - I was wondering how I can bring this functionality to Azure Functions, so that the script can be called via an API. As Pytho…

Best way to add python scripting into QT application?

I have a QT 4.6 application (C++ language) and i need to add python scripting to it on windows platform. Unfortunately, i never embed python before, and it seems to be a lot of different ways to do so.…

Unexpected behavior of python builtin str function

I am running into an issue with subtyping the str class because of the str.__call__ behavior I apparently do not understand. This is best illustrated by the simplified code below.class S(str):def __ini…

How to set cookies with GAE/Python for 1 month?

I need to implement the following:User input user id and pass We validate that on another server If they are correct, cookies with these details should be saved for one month Each time user uses my sit…

connection refused with Celery

I have a Django project on an Ubuntu EC2 node, which I have been using to set up an asynchronous using Celery. I am following How to list the queued items in celery? along with the docs, to experiment…